diff --git a/cara/model_scenarios.py b/cara/model_scenarios.py index a5d6d2a1..54fc199f 100644 --- a/cara/model_scenarios.py +++ b/cara/model_scenarios.py @@ -1,3 +1,4 @@ +from numpy.core.function_base import linspace from cara import models from cara.monte_carlo.data import activity_distributions from tqdm import tqdm @@ -8,6 +9,7 @@ from scipy.spatial import ConvexHull import pandas as pd import matplotlib.lines as mlines from matplotlib.patches import Rectangle +import matplotlib as mpl ######### Plot material ######### fig = plt.figure() @@ -128,6 +130,104 @@ def exposure_model_from_vl_talking(viral_loads): return er_means +def exposure_model_from_vl_talking_cn(viral_loads): + + n_lines = 5 + cns = np.linspace(0., .2, n_lines) + norm = mpl.colors.Normalize(vmin=cns.min(), vmax=cns.max()) + cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.jet) + cmap.set_array([]) + + for cn in tqdm(cns): + er_means = [] + er_medians = [] + lower_percentiles = [] + upper_percentiles = [] + for vl in viral_loads: + exposure_mc = mc.ExposureModel( + concentration_model=mc.ConcentrationModel( + room=models.Room(volume=100, humidity=0.5), + ventilation=models.AirChange( + active=models.SpecificInterval(((0, 24),)), + air_exch=0.25, + ), + infected=mc.InfectedPopulation( + number=1, + virus=models.Virus( + viral_load_in_sputum=10**vl, + infectious_dose=50., + ), + presence=mc.SpecificInterval(((0, 2),)), + mask=models.Mask.types["No mask"], + activity=activity_distributions['Seated'], + expiration=models.Expiration.types['Talking'], + ), + ), + exposed=mc.Population( + number=14, + presence=mc.SpecificInterval(((0, 2),)), + activity=models.Activity.types['Seated'], + mask=models.Mask.types["No mask"], + ), + ) + exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) + # divide by 4 to have in 15min (quarter of an hour) + emission_rate = exposure_model.concentration_model.infected.emission_rate_when_present(cn)/4 + er_means.append(np.mean(emission_rate)) + er_medians.append(np.median(emission_rate)) + lower_percentiles.append(np.quantile(emission_rate, 0.01)) + upper_percentiles.append(np.quantile(emission_rate, 0.99)) + + # divide by 4 to have in 15min (quarter of an hour) + coleman_etal_er_talking_2 = [x/4 for x in coleman_etal_er_talking] + + + ax.plot(viral_loads, er_means, color=cmap.to_rgba(cn)) + ax.fill_between(viral_loads, lower_percentiles, + upper_percentiles, alpha=0.2, color=cmap.to_rgba(cn)) + + fig.colorbar(cmap, ticks=cns) + ax.set_yscale('log') + + ############# Coleman ############# + plt.scatter(coleman_etal_vl_talking, coleman_etal_er_talking_2, marker='x') + x_hull, y_hull = get_enclosure_points( + coleman_etal_vl_talking, coleman_etal_er_talking_2) + # plot shape + plt.fill(x_hull, y_hull, '--', c='orange', alpha=0.2) + + ############# Markers ############# + markers = [5, 'd', 4] + + ############ Legend ############ + result_from_model = mlines.Line2D( + [], [], color='blue', marker='_', linestyle='None') + coleman = mlines.Line2D([], [], color='orange', + marker='x', linestyle='None') + + title_proxy = Rectangle((0, 0), 0, 0, color='w') + titles = ["$\\bf{CARA \, \\it{(SARS-CoV-2)}:}$", "$\\bf{Coleman \, et \, al. \, \\it{(SARS-CoV-2)}:}$"] + leg = plt.legend([title_proxy, result_from_model, title_proxy, coleman], + [titles[0], "Result from model", titles[1], "Dataset"]) + + # Move titles to the left + for item, label in zip(leg.legendHandles, leg.texts): + if label._text in titles: + width = item.get_window_extent(fig.canvas.get_renderer()).width + label.set_ha('left') + label.set_position((-3*width, 0)) + + ############ Plot ############ + plt.title('Exhaled virions while talking for 15min', + fontsize=16, fontweight="bold") + plt.ylabel( + 'Aerosol viral load, $\mathrm{vl_{out}}$\n(RNA copies)', fontsize=14) + plt.xticks(ticks=[i for i in range(2, 13)], labels=[ + '$10^{' + str(i) + '}$' for i in range(2, 13)]) + plt.xlabel('NP viral load, $\mathrm{vl_{in}}$\n(RNA copies)', fontsize=14) + plt.show() + + return er_means def exposure_model_from_vl_breathing(viral_loads): diff --git a/cara/models.py b/cara/models.py index 240bac76..35930a16 100644 --- a/cara/models.py +++ b/cara/models.py @@ -555,7 +555,7 @@ class Expiration(_ExpirationBase): BLO_factors: typing.Tuple[float, float, float] @cached() - def aerosols(self, mask: Mask): + def aerosols(self, mask: Mask, cn: float): """ Result is in mL.cm^-3 """ def volume(d): return (np.pi * d**3) / 6. @@ -565,9 +565,9 @@ class Expiration(_ExpirationBase): return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.262364)) * np.exp(-1 * (np.log(d) - 0.989541) ** 2 / (2 * 0.262364 ** 2))) - def _Lmode(d: float) -> float: + def _Lmode(d: float, cn: float) -> float: # L-mode (see ref. above). - return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.506818)) * + return ( (1 / d) * (cn / (np.sqrt(2 * np.pi) * 0.506818)) * np.exp(-1 * (np.log(d) - 1.38629) ** 2 / (2 * 0.506818 ** 2))) def _Omode(d: float) -> float: @@ -577,7 +577,7 @@ class Expiration(_ExpirationBase): def integrand(d: float) -> float: return (self.BLO_factors[0] * _Bmode(d) + - self.BLO_factors[1] * _Lmode(d) + + self.BLO_factors[1] * _Lmode(d, cn) + self.BLO_factors[2] * _Omode(d) ) * volume(d) * (1 - mask.exhale_efficiency(d)) @@ -670,7 +670,7 @@ class InfectedPopulation(Population): #: The type of expiration that is being emitted whilst doing the activity. expiration: _ExpirationBase - def emission_rate_when_present(self) -> _VectorisedFloat: + def emission_rate_when_present(self, cn: float) -> _VectorisedFloat: """ The emission rate if the infected population is present. @@ -680,7 +680,7 @@ class InfectedPopulation(Population): # Emission Rate (virions / h) # Note on units: exhalation rate is in m^3/h, aerosols in mL/cm^3 # and viral load in virus/mL -> 1e6 conversion factor - aerosols = self.expiration.aerosols(self.mask) + aerosols = self.expiration.aerosols(self.mask, cn) ER = (self.virus.viral_load_in_sputum * self.activity.exhalation_rate * diff --git a/cara/plot_output.py b/cara/plot_output.py index e7f58869..755ac12b 100644 --- a/cara/plot_output.py +++ b/cara/plot_output.py @@ -1,11 +1,12 @@ -from cara.model_scenarios import exposure_model_from_vl_breathing, exposure_model_from_vl_talking +from cara.model_scenarios import exposure_model_from_vl_breathing, exposure_model_from_vl_talking, exposure_model_from_vl_talking_cn import numpy as np import csv viral_loads = np.linspace(2, 12, 600) #er_means = exposure_model_from_vl_talking(viral_loads) -er_means = exposure_model_from_vl_breathing(viral_loads) +#er_means = exposure_model_from_vl_breathing(viral_loads) +er_means = exposure_model_from_vl_talking_cn(viral_loads) with open('data.csv', 'w', newline='') as csvfile: fieldnames = ['viral load', 'emission rate']