• Modeling the role of respiratory droplets in Covid-19 type pandemics
• 

For the modeling purpose, the exhaled droplets are assumed to evaporate in a quiescent environment at a fixed ambient temperature and relative humidity. In reality, during coughing, talking, or sneezing, the droplets are exhaled in a turbulent jet/puff.¹⁶16. L. Bourouiba, E. Dehandschoewercker, and J. W. Bush, “Violent expiratory events: On coughing and sneezing,” J. Fluid Mech. 745, 537–563 (2014). https://doi.org/10.1017/jfm.2014.88 However, as shown in Eqs. (10) and (11), the puff rapidly decelerates due to entrainment and lack of sustained momentum source, rendering the average V_D,P to be less than 1% of the initial velocity. Furthermore, since the Prandtl number, defined as ratio of kinematic viscosity and thermal diffusivity, is approximately unity (Pr = ν/α ≈ 0.71) for air, we can safely assume that the temperature and relative humidity that the droplets in the puff experience are on average very close to that of the ambient. At the initial stages, the puff will indeed be slightly affected by buoyancy, which will influence droplet cooling and evaporation dynamics. Quantifying these effects accurately, merit separate studies, see for e.g., Ref. 3232. R. Narasimha, S. S. Diwan, S. Duvvuri, K. R. Sreenivas, and G. S. Bhat, “Laboratory simulations show diabatic heating drives cumulus-cloud evolution and entrainment,” Proc. Natl. Acad. Sci. U. S. A. 108, 16164–16169 (2011). https://doi.org/10.1073/pnas.1112281108 for buoyant clouds. In a higher dimensional model, these could be incorporated. Nonetheless, the evaporation rate of the droplet is driven by the transport of water vapor from the droplet surface to the ambient far field. Assuming the quasi-steady state condition, the evaporation mass flux can be written as

$$\begin{array}{cc}\hfill {\dot{m}}_1& =-4\pi {\rho }_v{D}_v{R}_s\mathrm{l}\mathrm{o}\mathrm{g}\left(1+B_M\right),\hfill \\ \hfill {\dot{m}}_1& =-4\pi {\rho }_v{\alpha }_g{R}_s\mathrm{l}\mathrm{o}\mathrm{g}\left(1+B_T\right).\hfill \end{array}$$

(15)

Here, ṁ₁ is the rate of change of the droplet water mass due to evaporation, R_s is the instantaneous droplet radius, ρ_v is the density of water vapor, D_v is the binary diffusivity of water vapor in air, and α_g is the thermal diffusivity of surrounding air. B_M = (Y_1,s − Y_1,∞)/(1 − Y_1,s) and B_T = C_p,l(T_s − T_∞)/h_fg are the Spalding mass transfer and heat transfer numbers, respectively. Here, Y₁ is the mass fraction of water vapor, while subscripts s and ∞ denote the location at the droplet surface and at the far field, respectively. The numerical subscripts 1, 2, and 3 will denote water, air, and salt, respectively. C_p,l and h_fg are the specific heat and specific latent heat of vaporization of the droplet liquid. For the pure water droplet, the vapor at the droplet surface can be assumed to be at the saturated state. However, as indicated earlier, the exhaled droplets during talking, coughing, or sneezing are not necessarily pure water; rather, they contain plethora of dissolved substances.⁵5. X. Xie, Y. Li, A. Chwang, P. Ho, and W. Seto, “How far droplets can move in indoor environments—Revisiting the wells evaporation-falling curve,” Indoor air 17, 211–225 (2007). https://doi.org/10.1111/j.1600-0668.2007.00469.x The existence of these dissolved non-volatile substances, henceforth denoted as solute, significantly affects the evaporation of these droplets by suppressing the vapor pressure at the droplet surface. The modified vapor pressure at the droplet surface for binary solution can be expressed by Raoult’s Law, P_vap(T_s, χ_1,s) = χ_1,sP_sat(T_s), where χ_1,s is the mole-fraction of the evaporating solvent (here water) at the droplet surface in the liquid phase¹⁴14. W. A. Sirignano, Fluid Dynamics and Transport of Droplet and Sprays (Cambridge University Press, 2010). and χ_1,s = 1 − χ_3,s. The far field vapor concentration, on the other hand, is related to the relative humidity of the ambient. Considering the effects of Raoult’s law and relative humidity, the vapor concentrations at the droplet surface and at the far field can be expressed as

$$\begin{array}{cc}\hfill Y_{1,s}& =\frac{P_{vap}\left(T_s,{\chi }_{1,s}\right)M_1}{P_{vap}\left(T_s,{\chi }_{1,s}\right)M_1+\left(1-P_{vap}\left(T_s,{\chi }_{1,s}\right)\right)M_2},\hfill \\ \hfill Y_{1,\infty }& =\frac{\left(RH\right)P_{sat}\left(T_{\infty }\right)M_1}{\left(RH\right)P_{sat}\left(T_{\infty }\right)M_1+\left(1-\left(RH\right)P_{sat}\left(T_{\infty }\right)\right)M_2}.\hfill \end{array}$$

(16)

$M_1$ and $M_2$ denote the molecular weights of water and air, respectively. For evaporation, the droplet requires latent heat, which is provided by the droplet’s internal energy and surrounding ambient. It has been verified that the thermal gradient in the liquid phase is rather small. Therefore, neglecting the internal thermal gradients, the energy balance is given by

$$m{C}_{p,l}\frac{\partial T_s}{\partial t}=-k_g{A}_s\frac{\partial T}{\partial r}{|}_s+{\dot{m}}_1{h}_{fg}-{\dot{m}}_1{e}_l,$$

(17)

where T_s is instantaneous droplet temperature, $m=\left(4/3\right)\pi {\rho }_l{R}_s^3$ and $A_s=4\pi R_s^2$ are the instantaneous mass and surface area of the droplet, ρ_l and e_l are the density and specific internal energy of the binary mixture of salt (if present) and water, and k_g is the conductivity of gas surrounding the droplet. $\frac{\partial T}{\partial r}{|}_s$ is the thermal gradient at the droplet surface and can be approximated as (T_s − T_∞)/R_s, which is identical to convective heat transfer for a sphere with a Nusselt number of 2. As such, including aerodynamic effects, the Nusselt number is given by $Nu=2+0.6R{e}_p^{0.5}P{r}^{0.3}$. The droplet Reynolds number, Re_p, was observed to be mostly less than 0.1, and as such, the aerodynamic enhancement of the Nusselt number, i.e., the second term in the right-hand side, is ignored.

Evaporative loss of water leads to an increase in the salt concentration in the droplet with time. As shown before, P_vap(T_s, χ_1,s) is a function of the salt concentration in the droplet, which thus must be modeled using the species balance equation, as shown in the following equation:

$$\frac{dm{Y}_3}{dt}+{\dot{m}}_{3,out}=0.$$

(18)

Here, Y₃ is the dissolved solute (salt) mass fraction in the droplet. ṁ_3,out, which represents the rate at which solute (salt) mass leaves the solution due to crystallization, is modeled below. Clearly, Eq. (18) shows that as water leaves the droplet, Y₃ increases. When Y₃ is sufficiently large such that the supersaturation ratio S = Y₃/Y_3,c exceeds unity, crystallization begins. Here, we use Y_3,c = 0.393 based on the efflorescent concentration of 648 g/l reported for NaCl–water droplets in Ref. 3333. F. Gregson, J. Robinson, R. Miles, C. Royall, and J. Reid, “Drying kinetics of salt solution droplets: Water evaporation rates and crystallization,” J. Phys. Chem. B 123, 266–276 (2018). https://doi.org/10.1021/acs.jpcb.8b09584. The growth rate of the crystal could be modeled using a simplified rate equation from^34,3534. A. Naillon, P. Duru, M. Marcoux, and M. Prat, “Evaporation with sodium chloride crystallization in a capillary tube,” J. Cryst. Growth 422, 52–61 (2015). https://doi.org/10.1016/j.jcrysgro.2015.04.01035. H. Derluyn, “Salt transport and crystallization in porous limestone: Neutron-x-ray imaging and poromechanical modeling,” Ph.D. thesis, ETH Zurich, 2012.

$$\frac{dl}{dt}={\left(S-1\right)}^{g_{cr}}{C}_{cr}{e}^{-E_a/R{T}_s}.$$

(19)

Here, l is the crystal length. Following Ref. 3535. H. Derluyn, “Salt transport and crystallization in porous limestone: Neutron-x-ray imaging and poromechanical modeling,” Ph.D. thesis, ETH Zurich, 2012., for NaCl, we find the constant C_cr = 1.14 × 10⁴ m/s, the activation energy E_a = 58 180 J/mol, and constant g_cr = 1. Using this, the rate of change of the crystal mass, which equals ṁ_3,out, is given by³⁵35. H. Derluyn, “Salt transport and crystallization in porous limestone: Neutron-x-ray imaging and poromechanical modeling,” Ph.D. thesis, ETH Zurich, 2012.

$${\dot{m}}_{3,out}=\frac{d{m}_{3,crystal}}{dt}=6{\rho }_s{\left(2l\right)}^2\frac{dl}{dt}.$$

(20)

We note that while crystallization process could involve complex kinetics of solute, solvent, and ions; a well-studied³⁵35. H. Derluyn, “Salt transport and crystallization in porous limestone: Neutron-x-ray imaging and poromechanical modeling,” Ph.D. thesis, ETH Zurich, 2012. single-step crystallization kinetics has been used here for tractability. It will be shown that this model is able to predict the experimentally studied droplet lifetime reasonably well.

The governing equations [Eqs. (15)–(20)] manifest that several physical mechanisms are coupled during the evaporation process. For T_s,0 > T_∞, the droplet undergoes rapid cooling from its initial value. The droplet temperature, however, should eventually reach a steady state limit (wet bulb). This limit is such that the droplet surface temperature will be lower than the ambient, implying a positive temperature gradient or heat input. The heat subsequently transferred from the ambient to the droplet surface after attaining the wet-bulb limit is used completely for evaporating the drop without any change in sensible enthalpy. For a droplet with pure water, i.e., no dissolved non-volatile content, the mole-fraction of the solvent at the surface remains constant at 100%, and at the limit of steady state, the droplet evaporation can be written in terms of the well-known D² law,^14,1814. W. A. Sirignano, Fluid Dynamics and Transport of Droplet and Sprays (Cambridge University Press, 2010).18. C. K. Law, Combustion Physics (Cambridge University Press, 2006).

$$D_s^2\left(t\right)=D_{s,0}^2-K_{m}t,$$

(21)

where

$$K_m=8\left({\rho }_v/{\rho }_l\right)D_{v}ln\left(1+B_M\right).$$

(22)

However, for a droplet with the binary solution, the evaporation becomes strongly dependent on the solvent (or solute) mole-fraction, which reduces (or increases) with evaporative mass loss. The transient analysis, thus, becomes critically important in determining the evolution of the droplet surface temperature and instantaneous droplet size. During evaporation, the mole-fraction of the solute increases and attains a critical super-saturation limit, which triggers precipitation. The precipitation and accompanied crystallization dynamics, essentially, reduce the solute mass dissolved the in liquid phase, leading to a momentary decrease in its mole-fraction. This, in turn, increases the evaporation rate as mandated by Raoult’s law, which subsequently increases the solute concentration. These competing mechanisms control evaporation at the latter stages of the droplet lifetime. At a certain point, due to continuous evaporation, the liquid mass completely depletes and evaporation stops. The droplet after complete desiccation consists only of salt crystals, probably encapsulating the viruses and rendering them inactive. If the SARS-CoV-2 virus would remain active within the salt crystal, also known as droplet nucleus or aerosol, Covid-19 could spread by aerosol transmission in addition to that by droplets. In this paper, we focus on infection spread exclusive by respiratory droplets since the role of aerosols is not clear for transmission of Covid-19.

To validate the model, few targeted experiments were conducted to observe isolated levitated droplets evaporating in a fixed ambient condition. Particularly, the droplets with (1% w/w) NaCl solution vaporized to shrink to 30% of its initial diameter during the first stage of evaporation, as shown in Fig. 3. Hereafter, a plateau-like stage is approached due to increased solute accumulation near the droplet’s surface, which inhibits the diameter from shrinking rapidly. However, as shown in Fig. 3, shrinkage does occur (until D_s/D_s,0 ≈ 0.2) as the droplet undergoes a sol–gel transformation. The final shape of the precipitate is better observed from the micrographs presented in Fig. 3.

Figure 3 shows the final precipitate morphology for the desiccated droplets. The precipitates display a cuboid shaped crystalline formation, which is consistent with the structure of the NaCl crystal. The size and crystallite structure does show some variation, which could be linked with the initial size of the droplet. Only precipitates from larger droplets could be collected since smaller sized precipitates tend to de-stabilize and fly-off the levitator post-desiccation. While the precipitate from larger sized droplets tend to yield larger and less number of crystals, smaller droplets seem to degenerate into even smaller crystallites. However, this work does not investigate the dynamics of morphological changes of crystallization in levitated droplets.

Figure 3 also displays the comparison between results obtained from experiments and modeling. Experiments were performed with both pure water droplets as well as with droplets of 1% salt solutions. Experiments have been described in Sec. II. For the pure water cases shown in the left panel, simulation results follow the experiments rather closely. In the pure water case, classical D² law behavior could be observed. For the salt water droplets, a deviation from the D² law behavior occurs, and droplet evaporation is slowed. This is governed by Raoult’s law where the reduced vapor pressure P_vap(T_s, χ_1,s) on the droplet surface results from the increasing salt concentration with time. The evaporation rate approaches zero at about D_s/D_s,0 for 0.3 (experiments) and 0.25 (simulations), respectively. However, the salt concentration attained at this stage exceeds the supersaturation S ≥ 1 required for onset of crystallization. Thus, the salt crystallizes, reducing its concentration and increasing P_vap(T_s, χ_1,s) such that evaporation and water mass loss can proceed until nearly all the water has evaporated and only a piece of solid crystal as shown in Fig. 3 is left. It can be observed from Fig. 3 that in all cases, the final evaporation time is predicted within 15% of the experimental values. This suggests that despite the model being devoid of complexities (associated with inhomogeneities of temperature and solute mass fraction within the droplet and simple one step reaction to model the crystallization kinetics), it demonstrates reasonably good predictive capability. It is prudent to mention again that although we have done the analysis for the single isolated droplet, in reality, coughing or sneezing involves a whole gamut of droplet sizes in the form of a cloud.

Humans expel respiratory droplets while sneezing, coughing, or talking loudly. Such droplets have a size range from 5 μm to 2000 µm,³⁶36. J. Duguid, “The size and the duration of air-carriage of respiratory droplets and droplet-nuclei,” Epidemiol. Infect. 44, 471–479 (1946). https://doi.org/10.1017/s0022172400019288 while the dispersion could depend on the severity of the action. For example, while talking, an average human being will expel ∼600 droplets in a size range of 25 µm–50 µm, but this number goes upto ∼800 in the case of sneezing. For any given act (sneezing, coughing, or talking), the highest number of droplets fall in the range between 25 µm and 50 µm, while the smaller or larger droplets (than the above) are comparatively fewer in number. Nevertheless, the expelled volume of air contains a very small fraction of liquid droplets. To illustrate this point, the droplets are assumed to be in a uniform dispersion. The total volume of air expelled by a human being is estimated to be V_g ∼ 0.0005 m³. The total volume of liquid for a given mean droplet size is simply (NV_D,s,0), where N is the total number of droplets of size D_s,0 and V_D,s,0 is the volume of such a droplet. The total volume occupied by droplets of different sizes in a given act of coughing or sneezing is ${\sum }_i{N}_i{V}_{D_{s,0},i}$. Thus, the volume fraction (ϕ) of liquid droplets during a given act is ${\sum }_i{N}_i{V}_{D_{s,0},i}/V_g$. Using size distribution, reported in Ref. 3636. J. Duguid, “The size and the duration of air-carriage of respiratory droplets and droplet-nuclei,” Epidemiol. Infect. 44, 471–479 (1946). https://doi.org/10.1017/s0022172400019288, one can show that ϕ for sneezing, coughing, coughing with covered mouth, and talking loudly is 59 ppm, 549 ppm, 361 ppm, and 263 ppm, respectively. These numbers indicate that respiratory spray is rather sparse, implying that collective evaporation of droplet clusters may not be significant. This also justifies the modeling based on an isolated, contact free droplet.

Next, we set out to use this model to predict the droplet lifetime characteristics over a wide range of ambient conditions. The droplet evaporation time t_evap(D_s,0) is calculated from the analysis presented in Secs. V A and V B. Indeed, the droplet evaporation competes with gravitational settling.²2. W. Wells, “On air-borne infection: Study II. Droplets and droplet nuclei,” Am. J. Epidemiol. 20, 611–618 (1934). https://doi.org/10.1093/oxfordjournals.aje.a118097 The settling time t_settle(D_s,0) is calculated by accounting for the decreasing diameter using the equation for the Stokes settling velocity,

$$w=\left({\rho }_p-{\rho }_f\right)g{D}_s^2/18\mu .$$

(23)

The settling time is estimated as that time by which the droplet gets out of the radius from which breathable air is collected—already defined as σ_H/2 in Sec. III. Mathematically, t_settle is obtained by the following equation:

$${\int }_0^{t_{settle}}wdt={\sigma }_H/2.$$

(24)

Clearly, for any condition, while t_evap monotonically increases with D_s,0, t_settle monotonically decreases with D_s,0. In view of this, it is necessary to estimate the maximum time an exhaled droplet can remain within the collection volume without being evaporated or settled. Such a time can be estimated by defining a characteristic droplet lifetime τ, where

$$\tau =min\left\{t_{evap}|t_{evap}\ge t_{settle}\forall D_{s,0}\right\}.$$

(25)

τ is essentially the time where the two curves t_evap, t_evap as a function of D_s,0 intersect and represents the maximum time a liquid droplet of any size can exist before it is removed either by evaporation or gravity. For a given ambience specified by the ordered pair (T_∞, RH_∞), D_s,0 corresponding to τ can be defined as D_crit. Droplets with D_s,0 > D_crit settle due to gravity, while D_s,0 ≤ D_crit evaporate, earlier than the droplets with D_s,0 = D_crit. While droplets with D_s,0 ≠ D_crit can certainly transmit the disease, those with D_crit establishes the boundaries in terms of the lifetime, cloud diameter, and maximum distance traversed. D_crit is dependent on ambient conditions, i.e., temperature and relative humidity. The distribution of D_crit over a wide range of relevant ambient conditions is shown in Fig. 4(a). Interestingly, at high T_∞ and low RH_∞, where the evaporation rate is very fast, even a large droplet rapidly shrinks before it can settle. By the same argument at low T_∞ and high RH_∞, a relatively smaller droplet cannot evaporate quickly; therefore, t_evap = t_settle is attained for smaller droplet sizes. This explains why we observe large D_crit at high T_∞, low RH_∞, and small D_crit at low T_∞ and high RH_∞. The evaporation time of the droplet of diameter D_crit, which has been established as the characteristic lifetime τ of the droplet set, is shown in Fig. 4(b). Despite the different initial sizes, we find that τ is minimum at high T_∞ and low RH_∞ conditions, whereas it is maximum at low T_∞ and high RH_∞. At the same time, longer lifetime, i.e., large τ, allows the droplet cloud to travel a longer distance axially (X_p) and disperse radially (σ_D). Thus, X_P and σ_D are observed to be smaller (larger) for high (low) T_∞ and low (high) RH_∞ conditions, as shown in Figs. 4(c) and 4(d). This also shows how the minimum required “social distance” represented by X_p is not constant but depends on ambient conditions. Combining these results, it can be concluded that the size, lifetime, distance traveled, and the radial dispersion of the longest surviving droplet is not constant and is a strong function of ambient conditions. In particular, low temperature and high RH enhance the droplet lifetime significantly. Relative humidity strongly affects the droplet lifetime compared to temperature. An increase in droplet lifetime also implies that such droplets stay in the ambient for longer periods and hence travel longer distances as reflected by X_p. This implies that such droplets can lead to higher infection propensities. The critical size, droplet lifetime, distance traveled, or size of the droplet cloud for any practical condition are readily obtainable from Figs. 4(a)–4(d), respectively.

In Fig. 5, we look into the evolution of the normalized mass and temperature of the droplets for two cases named as case A and case B, also identified with black dots in Fig. 4. For case A, (T_∞, RH_∞) = (8, 55), while for case B, (T_∞, RH_∞) = (28, 77). The unit of T_∞ is °C and that of RH is %. These conditions have been specifically chosen to loosely represent the spring weather in North America and South East Asia, respectively, for the early days of the Covid-19 pandemic at both these locations. Figure 5 clearly explains why τ_A > τ_B. Indeed, higher RH at case B implies that the temporary hiatus in evaporation due to reduced vapor pressure is reached at a higher water mass load of the droplet than at the case A condition. In both cases, this occurs at about 3 s. However, the higher temperature in B results in faster crystallization kinetics due to the Arrhenius nature of the equation given by Eq. (19), which causes an eventual faster crystallization rate than at A. Figure 5 clearly shows that although the knee in the solvent depletion profiles are attained at the same time, it is the crystallization and simultaneous desiccation dynamics that governs the eventual difference in lifetime of the droplets at two conditions. It remains to be seen whether this result holds for a detailed crystallization reaction mechanism. The temperature evolution plots shown in the right panel of Fig. 5 also reveals how the droplet initially exhaled at 30 °C rapidly cools to the corresponding wet-bulb temperature to subsequently allow heat transfer into the droplet leading to evaporation. However, as the salt concentration reduces due to crystallization, the temperature rises subsequently above the corresponding wet-bulb limits.

With the droplet lifetime available over a wide range of conditions, the corresponding infection rate constant and eigenvalues given by Eqs. (5) and (6) could be evaluated. Just to recapitulate, τ determines the infection rate constant k₁ by Eq. (14). In turn, k₁ affects the exponents of the time dependent infection equation [Eq. (7)], the growth parameters—eigenvalues λ₁, λ₂ through Eq. (8). The contours of λ₁, λ₂, and k₁ as a function of T_∞ and RH_∞ are shown in Figs. 6(a)–6(c), respectively. The direct correspondence between τ and k₁ is immediately apparent upon comparing the respective Figs. 4(b) and 6(c). The infection rate constant is highest at low T_∞ and high RH_∞ where the droplet evaporation is slowed due to the slow mass loss rate and enhanced crystallization time. On the other hand, faster droplet evaporation leads to small infection rate constant values at high T_∞ and low RH_∞. The temperature and relative humidity dependency of the eigenvalues λ₁ and λ₂ are shown in Figs. 6(a) and 6(b), respectively. The direct correspondence of Figs. 6(a) and 6(b) with Fig. 6(c) and Fig. 4 are established through ${\lambda }_{\mathrm{1,2}}=-0.54\pm \sqrt{0.22+k_1}$ given by Eq. (8). It should, however, be noted that due to the inherent negative sign of λ₂, its influence on determining the growth rate of the infected population is rather limited. It is λ₁ that primarily drives the growth of the infected population as apparent from Eq. (7). From Fig. 6(a), we observe that for a fixed T_∞, λ₁ increases with RH_∞, while for a fixed RH_∞, λ₁ decreases with T_∞. Furthermore, we observe that the iso-λ₁ contour lines bend and converge at RH_∞ > 75%. This means that for RH_∞ > 75%, large λ₁ > 0.4 is expected over a wider range of temperatures between 5 °C < T_∞ < 20 °C. This is a manifestation of the greatly reduced evaporation potential—the difference between water vapor concentration on the droplet surface and in the ambient at high RH_∞ conditions. This is further reflected in the dramatic difference in the rate ratio—the ratio of the cumulative number of positive cases on a particular day to the cumulative number of positive cases seven days before: N_P/N_P,0 in Fig. 6(d). This figure has been arrived at by assuming the local population density to be 10 000 km⁻². Furthermore, we calculate t_c in Eq. (14) as t_c = 3600 × 24/N_exp. N_exp is the number of infecting expiratory events per person per day and is assumed to be 3 based on the coughing frequency of 0–16 in normal subjects.³⁷37. J. Hsu, R. Stone, R. Logan-Sinclair, M. Worsdell, C. Busst, and K. Chung, “Coughing frequency in patients with persistent cough: Assessment using a 24 hour ambulatory recorder,” Eur. Respir. J. 7, 1246–1253 (1994). https://doi.org/10.1183/09031936.94.07071246 We find that purely based on ambient conditions, implying all other factors have been held constant, the rate ratio can be different by an order of magnitude between (T_∞, RH_∞) = (5, 75) vs (35, 20). Practically, such a contrast might be less apparent in real data in which other important factors such as population density, social enforcement, travel patterns, and susceptible supply³⁸38. R. E. Baker, W. Yang, G. A. Vecchi, C. J. E. Metcalf, and B. T. Grenfell, “Susceptible supply limits the role of climate in the early SARS-CoV-2 pandemic,” Science eabc2535 (2020). https://doi.org/10.1126/science.abc2535 exert significant influence.