Merge branch 'feature/diameter_dependent_dynamics_and_deposition' into 'feature/scientific_model_update'

Model upgraded into fully diameter-dependent model

See merge request cara/cara!263

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.

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)