diff --git a/cara/models.py b/cara/models.py index 157eb256..67d4c62a 100644 --- a/cara/models.py +++ b/cara/models.py @@ -768,13 +768,29 @@ class ConcentrationModel: ventilation: _VentilationBase infected: _PopulationWithVirus + #: evaporation factor: the particles' diameter is multiplied by this + # factor as soon as they are in the air (but AFTER going out of the, + # mask, if any). + evaporation_factor: float = 0.3 + @property def virus(self): return self.infected.virus def infectious_virus_removal_rate(self, time: float) -> _VectorisedFloat: - # Particle deposition on the floor (value from CERN-OPEN-2021-04) - vg = 1.88e-4 # value corresponding to a diameter of 2.5 microns + # Equilibrium velocity of particle motion toward the floor + if (isinstance(self.infected, InfectedPopulation) + and isinstance(self.infected.expiration, Expiration)): + d = self.infected.expiration.diameter * self.evaporation_factor + vg = 1.88e-4 * (d / 2.5)**2 # see CERN-OPEN-2021-04 + # (velocity of 1.88e-4 corresponds to diameter of 2.5 microns) + else: + # model is not evaluated for specific values of aerosol + # diameters - we choose a single velocity value + # corresponding to that obtained with a diameter of 2.5 microns + # (geometric average of the breathing expiration distribution, + # taking evaporation into account, see CERN-OPEN-2021-04) + vg = 1.88e-4 # Height of the emission source to the floor - i.e. mouth/nose (m) h = 1.5 # Deposition rate (h^-1) @@ -952,8 +968,26 @@ class ExposureModel: #: The number of times the exposure event is repeated (default 1). repeats: int = 1 - #: The fraction of viruses actually deposited in the respiratory tract - fraction_deposited: _VectorisedFloat = 0.6 + def fraction_deposited(self): + """ + The fraction of viruses actually deposited in the respiratory tract. + From W. C. Hinds, New York, Wiley, 1999 (pp. 233 – 259). + """ + if not (isinstance(self.concentration_model.infected,InfectedPopulation) + and isinstance(self.concentration_model.infected.expiration,Expiration)): + # model is not evaluated for specific values of aerosol + # diameters - we choose a single "average" deposition factor, + # as in CERN-OPEN-2021-04. + fdep = 0.6 + else: + # deposition factor depends on aerosol particle diameter. + d = (self.concentration_model.infected.expiration.diameter + * self.concentration_model.evaporation_factor) + IFrac = 1 - 0.5 * (1 - (1 / (1 + (0.00076*(d**2.8))))) + fdep = IFrac * (0.0587 + + (0.911/(1 + np.exp(4.77 + 1.485 * np.log(d)))) + + (0.943/(1 + np.exp(0.508 - 2.58 * np.log(d))))) + return fdep def _normed_exposure(self) -> _VectorisedFloat: """ @@ -992,15 +1026,44 @@ class ExposureModel: return (np.array(self._normed_exposure()*aerosols).mean() * emission_rate/aerosols) + def _deposited_exposure(self) -> _VectorisedFloat: + """ + The number of virus per meter^3 deposited on the respiratory + tract. As in the exposure method, with sampled diameters the + aerosol volume and the deposited fraction, have to be put + in the deposited exposure before taking the mean, to obtain the + proper result (which corresponds to an integration on diameters). + """ + emission_rate = self.concentration_model.infected.emission_rate_when_present() + if (not isinstance(self.concentration_model.infected,InfectedPopulation) + or not isinstance(self.concentration_model.infected.expiration,Expiration) + or np.isscalar(self.concentration_model.infected.expiration.diameter) + ): + # in all these cases, there is no distribution of + # diameters that need to be integrated over + return (self._normed_exposure() * + self.fraction_deposited() * + emission_rate) + else: + # the mean of the diameter-dependent exposure (including + # aerosols volume, but NOT the other factors) has to be + # taken first (this is equivalent to integrating over the + # diameters) + mask = self.concentration_model.infected.mask + aerosols = self.concentration_model.infected.expiration.aerosols(mask) + return (np.array(self._normed_exposure() * aerosols * + self.fraction_deposited()).mean() * + emission_rate/aerosols) + def infection_probability(self) -> _VectorisedFloat: - exposure = self.exposure() + deposited_exposure = self._deposited_exposure() f_inf = self.concentration_model.infected.fraction_of_infectious_virus() inf_aero = ( self.exposed.activity.inhalation_rate * (1 - self.exposed.mask.inhale_efficiency()) * - exposure * self.fraction_deposited * f_inf + deposited_exposure * f_inf ) # oneoverln2 multiplied by ID_50 corresponds to ID_63. diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index f6b22054..bb3f283b 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -55,6 +55,7 @@ def shared_office_mc(): presence=concentration_mc.infected.presence, activity=models.Activity.types['Seated'], mask=concentration_mc.infected.mask, + host_immunity=0., ), ) @@ -97,6 +98,7 @@ def classroom_mc(): presence=concentration_mc.infected.presence, activity=models.Activity.types['Seated'], mask=concentration_mc.infected.mask, + host_immunity=0., ), ) @@ -129,6 +131,7 @@ def ski_cabin_mc(): presence=concentration_mc.infected.presence, activity=models.Activity.types['Moderate activity'], mask=concentration_mc.infected.mask, + host_immunity=0., ), ) @@ -163,6 +166,7 @@ def gym_mc(): presence=concentration_mc.infected.presence, activity=models.Activity.types['Heavy exercise'], mask=concentration_mc.infected.mask, + host_immunity=0., ), ) @@ -196,6 +200,7 @@ def waiting_room_mc(): presence=concentration_mc.infected.presence, activity=models.Activity.types['Seated'], mask=concentration_mc.infected.mask, + host_immunity=0., ), ) @@ -230,7 +235,8 @@ def skagit_chorale_mc(): presence=concentration_mc.infected.presence, activity=models.Activity.types['Moderate activity'], mask=concentration_mc.infected.mask, - ), + host_immunity=0., + ), ) @@ -304,6 +310,7 @@ def test_small_shared_office_Geneva(mask_type, month, expected_pi, presence=concentration_mc.infected.presence, activity=activity_distributions['Seated'], mask=concentration_mc.infected.mask, + host_immunity=0., ), ) exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE)