A comprehensive spatial-temporal infection model
AA COMPREHENSIVE SPATIAL - TEMPORAL INFECTION MODEL
A P
REPRINT
Harisankar Ramaswamy
Aerospace and Mechanical EngineeringViterbi School of EngineeringUniversity of Southern California [email protected]
Assad A. Oberai
Aerospace and Mechanical EngineeringViterbi School of EngineeringUniversity of Southern California [email protected]
Yannis C. Yortsos ∗ Mork Family Department of Chemical Engineering and Materials ScienceViterbi School of EngineeringUniversity of Southern California [email protected]
August 31, 2020 A BSTRACT
Motivated by analogies between the spreading of human-to-human infections and of chemical pro-cesses, we develop a comprehensive model that accounts both for infection (reaction) and for transport(mobility, advection and diffusion). In this analogy, the three different populations (susceptible, in-fected and recovered) of infection models correspond to three “chemical species”. Areal densities(people/area), rather than populations, emerge as the key variables, thus capturing the effect of spatialdensity, widely considered important, but ignored or under-represented in existing models. We deriveexpressions for the kinetics of the infection rates and for the important parameter R , that includeareal density and its spatial distribution. Coupled with mobility (through diffusion) the model allowsthe study of various effects.We first present results for a “batch reactor”, the chemical process equivalent of the SIR model.Because density makes R a decreasing function of the process extent, the infection curves aredifferent and smaller than for the standard SIR model, the difference increasing with R . We showthat the effect of the initial conditions (density of infected individuals) is limited to the onset of theepidemic, everything else being equal. The same invariance is obtained for infection imported intoinitially non-infected regions. We derive effective infection curves for a number of cases, includinga back-and-forth "commute" between regions of low (e.g. “home”) and high (e.g. “work”) R environments.We then consider spatially distributed systems. We show that diffusion leads to traveling waves,which in 1-D geometries (rectilinear or radial) propagate at a constant speed and with a constantshape, both of which are sole functions of R . The infection curves are slightly different than forthe batch problem, as diffusion mitigates the infection intensity, thus leading to an effective lower R . The dimensional wave speed is found to be proportional to the product of the square root of thediffusivity and of an increasing function of R , confirming the importance of restricting mobility inarresting the propagation of infection. We examine the interaction of infection waves under variousconditions and scenarios, and extend the wave propagation analysis to 2-D heterogeneous systems. K eywords Infection Models · Chemical Processes · Wave Propagation · COVID-19 ∗ Corresponding author a r X i v : . [ q - b i o . P E ] A ug PREPRINT - A
UGUST
31, 2020
Understanding the spread of infectious diseases, where infection is from human to human, is of significant interest,greatly relevant to epidemics, such as the one the world is experiencing today with COVID-19. A plethora ofepidemiologic models have been proposed, developed and tested [1, 2, 3]. These are based on three essential components:(i) Identifying different populations (typically, susceptible, infected, and recovered (or perished), and their extensions todemographic or health history subcategories); (ii) Describing ways by which the various populations come in proximitywith one another; and (iii) Postulating rates by which members from one population covert into another.Widely used is the celebrated SIR model [1, 4], which captures essential aspects of contagion using three populations:susceptible (denoted by S ), infected (denoted by I ) and recovered (including perished) (denoted by R ). The originalmodel, and many of its recent variants, account for items (i) and (iii) above, but not for (ii). More importantly, the SIRrepresentation is in terms of the total populations, rather than their areal densities (people/area), which are expected tobe the most important contagion variables. Such models also ignore spatial transfer, or the areal spreading of infections,that ultimately causes populations of different densities to interact with one another.Another category is based on computing individual trajectories of typical members of the various populations and onsubsequently postulating probabilities of collisions and infection [5]. These agent-based models implicitly account forspatial density effects through various mobility, propagation and collision rules. In particular, for specific individuals,whose trajectories and potential contacts are known, one can predict or retrace past infection paths. This forms the basisof contact tracing methods.A different and, we believe, more fundamental way to model the problem, is to follow a chemical reaction engineeringapproach, in which one can cast contagion spreading in terms of corresponding transport and reaction models. Thisapproach relies on the following analogies: populations map into chemical species; densities (specifically, arealdensities) into molecular concentrations; infection rates into chemical reaction rates (where mass action kinetics apply);and spatial transport into advective and diffusive (or dispersive) fluxes. Then, one can express the relevant population(species) balances using a partial differential equation description. It is the objective of this paper to create such amethodology for the description of the spread of an epidemic in general, with application to COVID-19, in particular.The work presented parallels related previous continuum models [6] but differs in a number of aspects to be explicitlydiscussed below.A key, fundamental underlying assumption is that we can homogenize population distributions by defining continuumvariables in terms of their areal averages. This allows one to postulate continuum conservation laws in terms of theirrate of change, transport and reaction. We recognize the limitations in this analogy between transport and interaction ofhuman populations on the one hand and molecular transport and reaction on the other. Given the relatively sparse arealdensities of human populations, compared to molecular densities, the law of large numbers may in fact not be in effect,in order to obtain well-defined averages (homogenization) [7], in which case the corresponding distributions might beakin to rarefied gas dynamics [8]. In addition, human movement is often influenced by behavioral drivers, rather thanrandom walks, and it is possible that, e.g. a diffusion or dispersion process would require a description different thanBrownian motion. These alternatives will be pointed out where appropriate. Despite such limitations, however, webelieve that the analogy pursued here is useful and instructive, as it allows one to obtain important new insights onfundamental spatial aspects of the problem.The commonly used SIR model arises from our formulation by integrating the relevant partial differential equationsover specific control areas (e.g. a building, a manufacturing plant, a school, a city, a state, a country, etc.), namely byconsidering a “batch reactor” model, where one assumes perfect mixing [9]. We examine its validity. Infection rates,depending crucially on spatial (areal) density, lead to spatial non-uniformities, while transport leads to infection waves.Both impact the stationary assumption of the SIR models. More importantly, we find that the key variable R is notconstant during the process, contrary to the common assumption. Instead, it is found to depend on the areal spatialdensity and to decrease as a function of the extent of the contagion. We consider the solution of spatially-dependentproblems by including diffusion in both 1-D (rectilinear or radial) and in 2-D geometries. We discuss the asymptoticemergence of traveling infection waves, the speed and shape of which are found to be independent of geometry, andalso correspond to a slightly smaller effective value of R , than for the case of the batch SIR model.The paper is organized as follows: We first proceed with constructing an equivalent chemical reaction and transportprocess. After recasting all conservation equations in dimensionless form, we derive the associated key parameter R .Then, we consider the solution of a number of specific problems, from a “batch reactor” model (the SIR equivalent), tothe propagation of infection in time and space (including both 1-D and 2-D domains). Infection waves are found todepend only on R , or its effective equivalent, and for batch or 1-D processes, to be independent of initial conditions.This lack of influence of initial conditions on the shape of the infection curve (in space and/or in time) is notable, as itsignals universal scaling properties, a property implicitly assumed in many SIR-type models.2 PREPRINT - A
UGUST
31, 2020
Consider the conservation of mass of the three quantities of interest, susceptible, infected, and recovered (includingperished), and associate to them three equivalent “chemical species”. The corresponding mass balances read ∂N i ∂T + ∇ · ( q N i ) = −∇ · ( D i ) + r i , i = I, S, R. (1)where N i is population density (number/area), q i (if any) is an advective velocity vector, D i is the diffusive (ordispersive) flux, and r i is the net reaction rate of species i , and converts populations into one another due to infectionand recovery, respectively. The advective velocity q i denotes a transport or mobility term. While we allow in (1) thethree velocities to be different, in the remainder we will take q i = q for all populations. This common velocity isassumed independent of N i or its gradients, although it can be a function of space and/or time. But we also note inadvance, that in the remainder we will not explicitly examine any effects of advection, included here only for the sakeof completeness.Next, define N S + N I + N R = ρ, (2)and D S + D I + D R = D (3)where ρ is the total population density (total number/area) and D is the total diffusive (or dispersive flux), to obtain ∂ρ∂T + ∇ · ( q ρ ) = −∇ · ( D ) (4)as the equation governing the evolution of total density. How to represent the diffusive (or dispersive) fluxes requires some discussion. For a typical Fick’s law type of diffusion[10], one may take D i = − D ∇ N i , (5)where we defined a constant diffusion (or dispersion) coefficient D , and assumed that all species are indistinguishableas far as their physical properties is concerned. Then, Equation (4) becomes ∂ρ∂T + ∇ · ( q ρ ) = ∇ · ( D ∇ ρ ) (6)suggesting that the total population density ρ , in addition to being advected, also diffuses in the direction of a negativespatial gradient. While possible and perhaps even likely in human dynamics, the suggestion that the overall density maydiffuse, is contrary to the common continuity equation for fluids, e.g. incompressible fluids [11]. If we consider analternative approach and take instead D i = − Dρ ∇ ( N i /ρ ) , (7)the more familiar form emerges ∂ρ∂T + ∇ · ( q ρ ) = 0 (8)Diffusive fluxes of the type (7) do in fact arise in random walks on non-uniform lattices, in the limit when the latticebecomes continuous [12, 13]. While either (5-6) or (7-8) can be taken to describe diffusive transport, in the examplesbelow we will only assume the more familiar continuity equation versions (7) and (8). This has the additional advantagethat in the absence of advection, the total density is time-independent and a stationary function of space, which is usefulfor the solution of many problems of interest, as the equations now become decoupled.It must be also noted that diffusion or dispersion can occur by different than a Brownian motion type random walk,e.g. by Levy flights (where the distribution of diffusion steps has a heavy tail, thus allowing for occasional, althoughrare, large steps) [14]. Modeling such motions can be handled by an integro-differential equation description throughfractional derivatives [15, 16]. While worth considering, such an approach will not be pursued here.3 PREPRINT - A
UGUST
31, 2020Diffusion coefficients can be estimated by considering the ratio of the square of the average radius of a random walkover an associated time interval. For instance, for an office type environment where random walks may have a meanradius of m , over an − hr period, one obtains D ≈ − m /s , which is about three orders of magnitude larger thanthat for molecular diffusion in gases. A different approach to evaluate dispersion at much larger scales is by inferencefrom the spatial spread of the epidemic. As we describe later, diffusion leads to infection waves, whose propagationspeed depends on the assumed diffusion coefficient. Values of the order of D ≈ m /s or higher are derived in orderto match the observed propagation velocities. Clearly, these very large values reflect an aggregate mix of transportactivities (including advection), here expressed only through a diffusion mechanism, which, therefore must be viewedonly as an approximation. Consider, next, contagion rates. Following our chemical reaction analogy, we postulate the following two chemicalreactions S + I → I (9) I → R (10)that convert the three species to one-another. The stoichiometric coefficient 2 in the RHS of (9) indicates that onemember of infected species I is produced as a result of its interaction with one member of susceptible species S . Inturn, species I , consumed in reaction (10), leads to species R , produced in one-to-one stoichiometry. Both reactions areirreversible.Applying mass action kinetics [17] in the reactions (9) and (10) provides expressions for the reaction rates in terms ofthe respective concentrations or densities r S = − KN S N I , r I = KN S N I − Λ N I , r R = Λ N I , (11)where we introduced the reaction rate constants K and Λ . Equations (11) are the same as those for the SIR model,except that here the rates are correctly expressed in terms of areal densities, rather than in terms of the total populations,which is inaccurately assumed in the typical models.We hasten to note that for a more fine-grained model, the reaction rates would also depend on additional demographicsand/or health conditions of the individual species. Such fine-graining is possible by further partitioning the populationsinto subgroups [18]. The corresponding description would then be in terms of an equivalent “multicomponent mixture”[19], where equations (1) are recast in terms of a larger species vector, with reactions (9)-(10) extended appropriately,possibly involving a product of a reaction matrix with the species vector. For simplicity, this generalization will not befurther considered here.Relevant to linear (“first-order”) reactions, Λ has dimensions of inverse time, a typical value being / days − [6, 20, 21]. The kinetic parameter expresses the rate at which infected individuals recover (or die), on average, and isintrinsic to the infected fraction. This is not the case for non-linear (e.g. “second-order”) reactions, like those involved inthe rate of generation of new infections (reaction (9)). As reflected in the kinetic parameter K , which has dimensions ofinverse (time × (number/area)), the infection kinetics will depend on the duration, method and type of human-to-humancontact, the protective gear (PPE) of susceptible and infected individuals, various biological and physiological variables,the ambient environmental conditions (room air conditioning), spatial distancing and other parameters. The mostcontrollable among these factors are the frequency and degree of of encounters (collisions) between individuals, as wellas the intensity of interaction.Possible approaches to estimating K include kinetic theory models (e.g. similar to Maxwell-Boltzmann models, wherethe kinetic parameter is inversely proportional to the molecular mean free path (average length between collisions),itself inversely proportional to density [22]). At least in some domain, K should be increasing with spatial density,something ignored in previous SIR models. Directly applying Maxwell-Boltzmann-type kinetic theories, however, isnot possible in the present context and must be refined: Human encounters are not elastic collisions, and typically lastover finite time intervals. More importantly, effects of spatial distancing, a recognized key to the kinetics of contagion,must be captured in any model. Indeed, it is by now widely accepted that infection rates are negligible for densitiesbelow a limiting value (corresponding to a separation of m or f t , as also supported by fluid mechanical models ofdroplet propagation, and recommended in health policy guidelines [23, 24, 25]).We incorporate such aspects by postulating the following dependence K ( ρ ) = (cid:26) , ρ < ρ K F ( ρ − ρ ρ − ρ ) , ρ ≥ ρ (12)4 PREPRINT - A
UGUST
31, 2020where F ( x ) is an increasing function of x , F (0) = 0 , F (1) ≈ . The value ρ ≈ . m − represents the density belowwhich the reaction rate is negligible, while the upper limit ρ ≈ m − corresponds to the maximum “packing” density.In (12) we have separated spatial distancing effects, included through ρ , from factors, such as biological, environmental,facial covering, etc., which enter through K only (with K decreasing substantially as facial covering is applied).Equation (12) captures both spatial distancing and biological and environmental effects. Upon more careful inspection,however, we realize that the density dependence must be further modified to not include the recovered fractioncomponent, since that fraction does not affect contagion (assuming immunity in all recovered individuals). Therefore,we must modify the way density enters in (12) by replacing it with ρ (1 − r ) (where we defined r = Rρ , see also below).The new expression (replacing K ( ρ ) with K ( ρ (1 − r )) ) reads K = K ( ρ, r ) = (cid:26) , ρ (1 − r ) < ρ K F ( (1 − r ) ρ − ρ ρ − ρ ) , ρ (1 − r ) ≥ ρ (13)When the rate of infection is relatively low (and r (cid:28) ), the correction is negligible. For strong infection rates, on theother hand, it can be quite significant, as will be demonstrated later in the paper.Some additional remarks are pertinent: The above assume that an infected individual can infect a susceptible one atthe same constant rate. This is true either for asymptomatic, infected individuals, or for those who do not exhibitsymptoms until a few days following infection. It is not true, however, when infected individuals are isolated.Nonetheless, the previous formalism still applies: Assume that the infected species N I is further subdivided into onecategory containing asymptomatic individuals (denoted by A , with corresponding density N A ) and another containingquarantined individuals (denoted by Q , with corresponding density N Q ). The associated infection reaction rate is now KN S N A . Based further on the reasonable assumption that the percentage of those infected, but are asymptomaticor with mild symptoms, is a fixed fraction (e.g. a ) of the total fraction of the infected, the infection rate becomes KaN S N I . This is of the same dependence as before, hence the previous holds, except that the kinetic constant is nowmodified by parameter a . The same reasoning (with an additional parameter included) holds when infected individualsare contagious, but not identified as such, until sometime after infection. On the other hand, our approach cannot aseasily account for correlations between infected and susceptible individuals, for example when susceptible individualshave increased contact, hence higher probability of infection, with specific infected individuals related to them, e.g. byfamily, work or other proximity relations.A final remark relates to the practice of reporting area-wide averages (e.g. for states or countries). Given that almost allareas will never on average reach the minimum density required for infection (e.g. . m − ), area-wide averages oversubstantially heterogeneous density distributions are not very informative. Rather, distinguishing high-density areas(e.g. urban places, stadiums, schools, retirement homes, etc.) from low-density ones (e.g. farms, rural) is essential, withthe reporting of much more fine-grained area statistics being much more informative. Connecting the transmission ofinfection between areas of different density is possible using the present formalism, which includes spatial transportand/or diffusion, and where K is space-dependent, as further discussed below. We next proceed with a dimensionless notation. Denote dimensionless variables by lower case symbols and normalizespecies densities by ρ , time by / Λ , space by a characteristic external length scale l , K by K , and velocities by acharacteristic velocity q . Using equation (7) for diffusion we obtain ∂s∂t + Da ∇ · ( v s ) − C ∇ (ln ρ ) · ∇ s = ∇ · ( C ∇ s ) − R ( ρ, r ) si (14) ∂i∂t + Da ∇ · ( v i ) − C ∇ (ln ρ ) · ∇ i = ∇ · ( C ∇ s ) + R ( ρ, r ) si − i (15) ∂r∂t + Da ∇ · ( v r ) − C ∇ (ln ρ ) · ∇ r = ∇ · ( C ∇ r ) + i (16) ∂ρ∂t + Da ∇ · ( v ρ ) = 0 (17)Here, we defined the dimensionless Damkohler number Da = ql Λ , the dimensionless diffusion number C = D Λ l = φ − ,and the rescaled velocity v . Note that in the chemical reaction engineering literature, φ is known as the Thiele modulus[9]. From the dimensionless formulation arises a most important parameter, and one most commonly associated withthe SIR model, namely R ( ρ, r ) = K ρ Λ κ ( ρ, r ) (18)5 PREPRINT - A
UGUST
31, 2020where κ ( ρ, r ) = (cid:26) , (1 − r ) ρ < ρ F ( (1 − r ) ρ − ρ ρ − ρ ) , (1 − r ) ρ ≥ ρ (19)Equation (18) shows that R ( ρ, r ) is not constant, as typically assumed, but also depends both on the areal densityand on the extent of the process r . For example, the ratio of the final value R ( ρ, r ∞ ) to its initial R ( ρ, can beapproximated (assuming a linear function for F ) by R ( ρ, r ∞ ) R ( ρ,
0) = (1 − r ∞ ) ρ − ρ ρ − ρ ≈ − r ∞ (20)assuming ρ /ρ (cid:28) . For consistency, in the remainder, we will express the R dependence of various solutions to theproblem, e.g. r ∞ , through its value at the onset of the process, R ( ρ, . The finding that R ( ρ, r ) is not constant isconsistent with spatial distancing health policy guidelines. To our knowledge this property has not been noted before,even though it has important implications to the prediction of infection results, as shown below.For completeness, we also remark that had we used the different version for diffusion (5), instead, then equations (16),(15) and (16) would still remain in effect, subject to multiplying by a factor of the third term on the LHS of eachequation. Consider, now, the application of the previous formulation to special cases, starting first with the zero-dimensional(“batch reactor”) perfect-mixing problem.
Assume a “batch reactor” (closed model, with no input or output) and conditions of perfect mixing, in which all spatialdensities only depend on time. By integrating the differential equations (14)-(16) and using (˙) to denote time derivatives,we obtain the following system of ordinary differential equations ˙ s ( t ) = − R ( ρ, r ) si (21) ˙ i ( t ) = R ( ρ, r ) si − i (22) ˙ r ( t ) = i (23)subject to the closure s + i + r = 1 (24)and the initial conditions i (0) = i , s (0) ≡ s = 1 − i , r (0) = 0 . (25)Even though the overall density is time-independent, in such an SIR-like model, spatial effects enter through the effectof the extent of contagion on R (equations (21), (22)).Equations (21)-(23) produce non-trivial results when an initial, even infinitesimally small, seed of infected individuals( i ) is present. An analytical solution is possible. Substitute (23) into (21) and integrate to give: s = s exp( − (cid:90) r R ( ρ, r (cid:48) ) dr (cid:48) ) (26)thus, i = 1 − r − s exp( − (cid:90) r R ( ρ, r (cid:48) ) dr (cid:48) ) (27)Further substitution into (23) gives ˙ r ( t ) = 1 − r − s exp( − (cid:90) r R ( ρ, r (cid:48) ) dr (cid:48) ) (28)which can be integrated to the final solution t = (cid:90) r du − u − s exp( − (cid:82) u R ( ρ, r (cid:48) ) dr (cid:48) ) . (29)6 PREPRINT - A
UGUST
31, 2020We remark, in passing, that by expanding the exponential in (28) in a Taylor series, assuming constant R and keepingthe first three terms in the expansion, leads to a Riccati equation ˙ r ( t ) = 1 − r (1 − s R ( ρ, − s R ( ρ, r (30)accurate for small R . Such a model has been used as another alternative to the SIR description [26, 27, 28]. R Results from the solution of the above equations are shown in Figure 1. We first remark that from Equation (22), R ( ρ,
0) = 1 is the boundary demarcating two regions, where an initial infection either decays ( R ( ρ, < ) orgrows ( R ( ρ, > ). We will focus on the latter case ( R ( ρ, > ). Plotted in Figure 1 is the time variation of thethree different populations for R ( ρ,
0) = 2 . , as well as of the curves obtained under the SIR assumption of constant R (taken throughout the process to be equal to R ( ρ,
0) = 2 . ). The infection curves are of a similar shape, butsignificantly larger when the constant R value is used.Figure 1: Variation of susceptible ( s ), infected ( i ) and recovered ( r ) fractions as a function of non-dimensional timefor the case when R varies following equation (19) (solid lines) with R ( ρ,
0) = 2 . , and for the case of a constantvalue of R = R ( ρ,
0) = 2 . (dashed lines), which is the SIR assumption. Note the substantial difference between theresults obtained. The initial infected fraction is i (0) = 10 − .Of importance is the concept of herd immunity , denoted here by h = r ∞ , and defined as the asymptotic value of therecovered individuals, r ∞ . The latter is the solution of the equation − r ∞ − s exp( − (cid:90) r ∞ R ( ρ, r (cid:48) ) dr (cid:48) ) = 0 . (31)Results are shown in Figure 2. We note that h is an increasing function of R ( ρ, (with h → i as R ( ρ, → ).Also plotted in the same figure is the herd immunity calculated under the SIR assumption of a constant R = R ( ρ, .The corresponding values are significantly higher, even for relatively mild rates of infection. An effective constant R can be defined using the following average R ≡ (cid:82) r ∞ R ( ρ, r (cid:48) ) dr (cid:48) r ∞ (32)to relate an effective constant R to R ( ρ, . This relationship can be also inferred from Figure 2, which for exampleshows that to a constant value of R = 2 . corresponds a twice as large value of R ( ρ,
0) = 4 . . Plotted in Figure 2 isalso the maximum in the infected fraction, i max , which is an increasing function of R ( ρ, , as expected.A related quantity is the length of the epidemic, here approximated by t ∞ = (cid:90) r e du − u − s exp( − (cid:82) u R ( ρ, r (cid:48) ) dr (cid:48) ) . (33)7 PREPRINT - A
UGUST
31, 2020where r e = 0 . r ∞ , and r ∞ is the solution of (31). Results are shown in Figure 2. The dependence of t ∞ on R ( ρ, is monotonically decreasing, with smaller values of R ( ρ, (and lower infection levels) resulting into a longer duration,as long as 30 in dimensionless time (corresponding roughly to 15 months) or longer. A more typical duration, but withhigher infection rates, is about 12 dimensionless time units (6 months). The fact that the epidemic is higher at smallervalues of R is counter-intuitive, and calls for the need to educate the public that epidemics on the border of beingunder control will last longer.Figure 2: Variation of herd immunity ( r ∞ ), maximum infection fraction ( i max ) and the duration of a pandemic ( t ∞ ) asa function of R ( ρ, . Plotted also is the herd immunity ( r c, ∞ ) calculated assuming a constant R . The initial infectionfraction is i (0) = 0 . .Of relevance is the stability of the asymptotic state. Consider a small resurgence of infections after an asymptotic stateis reached. As long as condition R ( ρ, r ∞ )(1 − r ∞ ) < is satisfied (which is always the case), such fluctuations willdecay exponentially fast, as shown in equation (22) demonstrating that the final state is asymptotically stable. Suchwill not be the case, however, when the fluctuation is instead in R ( ρ, r ∞ ) , e.g. when a new value, e.g. R (cid:48) , such that R (cid:48) (1 − r ∞ ) > , sets in. For example, this could be the result of relaxations to spatial distancing, and of abandoningdue caution, after the asymptotic state is reached. Then, there will be an eruption of new cases, resulting into a newspreading of infections (a “second wave”), which will in turn reach a new asymptotic state, and a correspondingly new,higher herd immunity. Figure 3 demonstrates the development of such a second wave. The figure shows three differentregimes: An initial one with R ( ρ,
0) = 3 for which infection grows; a second regime following the imposition ofrestrictions at t = 1 (which leads to R = 0 . ), and which results in the "flattening of the curve" with a correspondingherd immunity of . ; and a regime in which relaxation of restrictions at t = 10 and a return to a higher R ( ρ,
0) = 3 ,leads to a second wave. The new epidemic lasts at least as long as the first one, and contributes to substantial newinfections (almost double the initial). This lack of structural stability of the asymptotic state is derives from the fact thatthe infection reaction (9) is autocatalytic, and leads, at least initially, to exponential rises. It illustrates the importance ofclosely adhering to a consistent, and as low as possible, value of R , for a desired herd immunity to be sustained. The autocatalytic nature of reaction (9) raises the additional question as to whether or not the infection curves dependon parameters other than R . Consider, first, the impact of the initial condition i . Figure 4 is a plot of the infectioncurve for different values of the initial condition for R ( ρ,
0) = 2 . . For the typical values considered, we observe thata decrease in the initial condition fraction leads to a shift in the infection curve i ( t ) to the right, but otherwise producesshapes that remain invariant. The shift is approximately 2 non-dimensional time units for each decrease in the initialcondition by a factor of 10.To explain these findings we consider the early-time solution of equation (22). At small times we obtain i ( t ) ≈ i o exp(( R ( ρ, − t ) = exp(( R ( ρ, − t − t )) , (34)8 PREPRINT - A
UGUST
31, 2020Figure 3: Variation of susceptible ( s ), infected ( i ) and recovered ( r ) fractions as a function of non-dimensional timewith R ( ρ,
0) = 3 , t ∈ (0 , , R ( ρ,
0) = 0 . , t ∈ (1 , , R ( ρ,
0) = 3 , t ∈ (10 , . The initial infected fraction is i (0) = 0 . .where t = − log( i ) ln(10) R ( ρ, − ≈ − . log( i ) R ( ρ, − . This confirms the existence of a time shift t of the same magnitude as inthe figure, consistent with the simulations. For the same reasons, the maximum infection fraction i max is independentof the initial condition i , assuming that the latter remains relatively small. The invariance in the shape of the infectioncurve and the fact that a decrease in the number of initial infections only acts to delay the onset of infection, aresignificant from a health policy perspective: containing the initial number of infections helps to provide a non-trivialcushion of time to contain the infection and educate the public to modify behavior by lowering R , thus ultimatelymitigating the intensity of the incoming infection. Absent such preparation or behavior modification, will negate anybeneficial effect of the lower number of initial infections: A corresponding contagion wave will emerge, ultimatelydependent only on R .Figure 4: Variation of infected ( i ) fraction as a function of non-dimensional time for different values of the initialinfected fraction i and R ( ρ,
0) = 2 . . 9 PREPRINT - A
UGUST
31, 2020
The same conclusions apply for “imported infection”. With this term, we refer to the case where a control area (e.g. acountry), originally without any infected individuals, receives a constant influx of infected and susceptible individualsover a finite (and small, reflecting a prompt public authority response) time interval τ , following which such influx isstopped (e.g. via a “flight ban”). By integrating equation (15) over spatial dimensions and assuming perfect mixing wefind to leading order ˙ i ( t ) = ( R ( ρ, − t + j i (35)valid for < t < τ . Here, we denoted the influx rate by j i and also assumed that the corresponding changes in densityare negligible, as long as τ is small. Solving equation (35) subject to the initial condition i (0) = 0 , we obtain, toleading-order (small τ ) i = j i t + · · · (36)It follows that as long as R ( ρ, τ (cid:28) j i , the impact of “imported infections” is simply to nucleate the process,essentially providing the equivalent initial condition, i (0) = i = j i τ , and the problem reverts to the previous. Weconclude that the net impact of a ban is to provide additional time (e.g. a month) for the necessary public healthmeasures to be implemented and to possibly lead to as small a value of R ( ρ, as possible, thereby limiting the spreadof the inevitable infection. In the absence of additional such prophylactic measures, a ban can only serve to delay theonset of contagion. The batch (SIR-like) model rests on the assumption that all profiles are spatially uniform. We explore the validity ofthis assumption in the presence of fluctuations or other non-uniformities, by considering the solution of a typical 1-Dproblem. We take ∂s∂t − C ∂ ln ρ∂x ∂s∂x = C ∂ s∂x − R ( ρ, r ) si (37) ∂i∂t − C ∂ ln ρ∂x ∂i∂x = C ∂ i∂x + R ( ρ, r ) si − i (38) ∂ρ∂t = 0 (39)subject to no-flux conditions at the two ends, ∂s∂x = ∂i∂x = ∂ρ∂x = 0 , at x = { , } . We consider the case where the initialfluctuations are variable, but the density is uniform, i ( x,
0) = i (1 + (cid:15)g ( x )) , s ( x,
0) = 1 − i ( x, , ρ ( x, t ) = ρ m , (40)or where the non-uniformity is in the density profile only i ( x,
0) = i , s ( x,
0) = 1 − i ( x, , ρ ( x, t ) = ρ m (1 + (cid:15)g ( x )) , (41)where (cid:15) (cid:28) . Expecting that diffusion ( C > ) will help to smooth spatial non-uniformities, we will focus in thissection on the solution in the absence of diffusion ( C = 0 ). Thus, any effects to arise will correspond to the spatialaveraging of an ensemble of batch problems. When C = 0 the relevant equations revert to (21)-(23). Consider a small perturbation expansion and denote means bysuperscript bar and fluctuations by superscript prime to obtain ∂ ¯ s∂t = − R ( ρ m )(¯ s ¯ i + s (cid:48) i (cid:48) ) (42) ∂ ¯ i∂t = + R ( ρ m )(¯ s ¯ i + s (cid:48) i (cid:48) ) − ¯ i (43)where we ignored the effect of the fluctuations on R . By further taking the representation s (cid:48) = (cid:15)i σ ( t ) g ( x ) and s (cid:48) = (cid:15)i η ( t ) g ( x ) , the fluctuations to leading order satisfy dσdt = − R ( ρ m )(¯ sη + ¯ iσ ) (44) dηdt = + R ( ρ m )(¯ sη + ¯ iσ ) − η (45)10 PREPRINT - A
UGUST
31, 2020with corresponding initial conditions σ (0) = − and η (0) = 1 .We first note that fluctuations contribute to the rate expression, Equations (42)-(43), through the mean of the product s (cid:48) i (cid:48) . This is negative and proportional to the square of the amplitude of the fluctuations ( i (cid:15)g ( x )) . A first conclusion,therefore, is that fluctuations lower the effective rate of the infection reaction (although only to order (cid:15) assuming thatthe fluctuations remain bounded). Adding diffusion will further reduce any such impact.To determine whether or not the fluctuations are bounded, we must find the eigenvalues ω of the matrix of the linearsystem (44)-(45), namely the solutions of ω + ( R ( ρ m )¯ i − R ( ρ m )¯ s + 1) + R ( ρ m )¯ i = 0 . (46)The equation has two negative eigenvalues as long as R ( ρ m )¯ i − R ( ρ m )¯ s + 1 > . (47)Equations (21)-(23) show that for some time before the infection fraction i reaches its maximum, and always after thattime, condition (47) is indeed valid. We conclude that the fluctuations will remain bounded, if not altogether decay,hence they have little effect on the average behavior. This finding is demonstrated in Figure 5, where we plot ensembleaverage responses. The curves obtained are practically the same either in the presence or the absence of fluctuations.Figure 5: Variation of the ensemble average of the infected ( i ) fraction as a function of non-dimensional time for R ( ρ,
0) = 2 . . The initial condition is equation (40) with i = 0 . , g ( x ) = sin( x ) , and (cid:15) = 0 . . Results are quite different when the spatial heterogeneity is larger, however. Here, the interpretation of the compositeaverage must be done with a careful understanding of the underlying heterogeneities.Consider, first, two regions with the same density, but with different initial conditions that differ by an order ofmagnitude (e.g. i = 10 − and i = 10 − ) i ( x,
0) = i H ( ξ − x ) + i H ( x − ξ ) (48)where < ξ < and H is the Heaviside step function. From the previous analysis we expect that the region withhigher initial infections will respond faster, the other response trailing by a time lag (e.g. see Figure 4 and Equation(34)). Accordingly, the composite behavior will be controlled initially by the region with the larger number of initialinfections, but at later times by the second region. Figure 6 shows that this is indeed the case.More interesting, perhaps, is the case where the area of interest consists again of two regions, e.g. "urban" in the interval < x < ξ and "suburban" in the interval ξ < x < We are interested in the infection curve, under the conditions ofa “commute” between the two regions: During the day, and for a certain time interval of duration λ (expressed as a11 PREPRINT - A
UGUST
31, 2020Figure 6: Variation of the infected ( i ) fraction as a function of non-dimensional time for the two different regions withdifferent initial infections and the composite curve ( R ( ρ,
0) = 2 . ).fraction of the 24-hr period), the density in the urban region, ρ u , is high, and in the suburban region it is negligible;while during the remaining part of the day and at night, the density in the urban region is zero, and in the suburban is ρ s .These two densities are related by a conservation equation ρ s = ρ u ξ − ξ . This setting is intended to model the commutebetween two regions (“home” and “work”, or “home” and “school”) where we expect significantly different R values.By further taking ξ (cid:28) , we will simply assume R = 0 for conditions at “home”.This problem can be solved based on the detailed transport and reaction equations derived earlier, that include advectionand diffusion. A simpler alternative, on the other hand, is to represent the problem as two "batch reactors", whose R values oscillate between the two values, R ,w ( ρ, r ) and , when the population is at “work” (or “school”) or at“home”, respectively. The results of this hypothetical “commute”, for an 8-hr work period ( λ = 1 / ), are shown inFigure 7. Superimposed are also the results corresponding to the value of the "work" parameter R ,w ( ρ, r ) appliedeverywhere (equivalently, by taking λ = 1 ). As expected, while commuting leads to a much lower effective value R ,eff ( ρ,
0) = 1 . , it does not suppress infection, for this example. As long as R ,eff ( ρ, > infection will occur,although by reducing exposure (decreasing λ ) the effect on decreasing the resulting effective rates is non-trivial. Indeed,we have found (not shown here for the sake of brevity) that for all value of λ examined, the corresponding effectiveis roughly equal to the arithmetic mean, R ,eff ( ρ, ≈ λR ,w . This finding suggests that there is a critical exposurevalue λ crit = R ,w , below which contagion is suppressed. These findings are new, and must also be interpreted withdue caution, and subject to the assumptions made. The preceding sections dealt with applications in zero-dimensional space (batch reactors or their ensembles). Inproblems that involve spatial dimensions, transport (via advection, diffusion or dispersion) plays two different roles:One is stabilizing, reducing the effect of small spatial non-uniformities, as discussed above. The other is in the oppositedirection, however, and helps spread the infection spatially. This section will only consider mobility effects due todiffusion.For a constant diffusion coefficient, in the absence of advection ( Da v = ), the relevant equations become ∂s∂t − C ∇ (ln ρ ) · ∇ s = C ∇ s − R ( ρ, r ) si (49) ∂i∂t − C ∇ (ln ρ ) · ∇ i = C ∇ i + R ( ρ, r ) si − i (50) ∂r∂t − C ∇ (ln ρ ) · ∇ r = C ∇ r + i (51)We will first consider problems in one spatial dimension. 12 PREPRINT - A
UGUST
31, 2020Figure 7: The infection curves for the case of a commute between two regions of high and low values of R ( R ,w ( ρ,
0) = 5 ) and exposure frequency of λ = 1 / (solid lines). Plotted also are the infection curves correspondingto a uniform R ,w ( ρ,
0) = 5 (equivalently to λ = 1 ) (dashed lines). The initial infection fraction is i (0) = 0 . . Assume 1-D rectilinear geometries, with spatially constant density. The initial infection conditions are non-zero in aspecified interval, − < x < , the rest of the domain being free of infections. We are interested in exploring howdiffusion leads to the spreading of the contagion from the infected region. Specifically, we will explore if travelingwaves develop at large times.To this effect, we introduce moving coordinate ξ = x − V t , where V is the wave velocity. Assuming that a steady-stateis reached in these coordinates (denoted by tilde) we obtain, − V ∂ ˜ s∂ξ = C ∂ ˜ s∂ξ − R ˜ s ˜ i (52) − V ∂ ˜ i∂ξ = C ∂ ˜ i∂ξ + R ˜ s ˜ i − ˜ i (53) − V ∂ ˜ r∂ξ = C ∂ ˜ r∂ξ + ˜ i (54)subject to no-flux conditions at the ends ∂ ˜ s∂ξ = ∂ ˜ i∂ξ = ∂ ˜ r∂ξ = 0 , at ξ = ±∞ . (55)We first observe that the invariance of (52)-(54) to the transformation ξ → − ξ , V → − V , suggests that there will betwo asymptotic waves, one moving to the right, with velocity V , and one to the left, with velocity − V . Let ξ u and ξ d betwo locations sufficiently upstream and and downstream, respectively, of the wave fronts. We then expect ˜ r ( ξ u ) = r V, ∞ , ˜ i ( ξ u ) = 0 , and ˜ r ( ξ d ) = 0 , ˜ i ( ξ d ) = 0 (56)where we have anticipated that r V, ∞ might not be identical to the r ∞ of the batch reactor result (equation (31)).The wave velocity can be determined by integrating (54) between ξ u and ξ d , V = 1 r V, ∞ (cid:90) ξ d ξ u ˜ idξ. (57)and, equivalently, V = 1 r V, ∞ (cid:90) ∞−∞ ˜ idξ. (58)13 PREPRINT - A
UGUST
31, 2020We find, therefore, that a wave solution exists, as long as the integral in (58) is non-zero, which is true if there isa contagion. This contagion wave spreads, driven by diffusion, with a velocity that expresses the intensity of thecontagion.
In the above, we can explicitly remove the C -dependence by introducing rescaled space coordinates and velocities, ξ = √ Cζ and V = W √ C . The equations remain the same with C set as C = 1 . Namely, the profiles, now, denoted bysuperscript hat, satisfy − W ∂ ˆ s∂ζ = ∂ ˆ s∂ζ − R ˆ s ˆ i (59) − W ∂ ˆ i∂ζ = ∂ ˆ i∂ζ + R ˆ s ˆ i − ˆ i (60) − W ∂ ˆ r∂ζ = ∂ ˆ r∂ζ + ˆ i (61)subject to no-flux conditions at the two ends. The velocity becomes W = 1 r V, ∞ (cid:90) ζ d ζ u ˆ idζ, (62)in dimensional form, V = √ D Λ r V, ∞ (cid:90) ζ d ζ u ˆ idζ. (63)and by extension W = 1 r V, ∞ (cid:90) ∞−∞ ˆ idζ, (64) V = √ D Λ r V, ∞ (cid:90) ∞−∞ ˆ idζ. (65)We note that although the initially infected region x ∈ ( a, b ) transforms into ζ ∈ ( aC − / , bC − / ) , thus containinga C -dependence, the asymptotic traveling wave is independent of initial conditions, and hence of C . This result issignificant: it shows that the limit C → is a singular one, which means that the batch reactor (SIR) problem is alsoa singular limit. The question then becomes, how different are the resulting curves in the presence of diffusion, andwhether or not one can define an equivalent, effective, value of R . These are explored below.Numerical results of the solution of (49)-(51) are shown in Figure 8, for a problem in which the initial infection regionis near the left boundary, and R ( ρ,
0) = 2 . . We observe that the infection profile evolves as a function of time, andreaches the asymptotic traveling wave after about t = 3 . The wavelike nature of the solution is clear if one plots theinfected fraction in space-time coordinates (Figure 9), where a ridge with a constant slope is clearly seen.Consider now, the answer to the question in how different are the infection profiles with diffusion included. Figure10 shows the relevant wave profiles at a fixed value of x (taken at x = − ), calculated using the full system ofequations, with diffusion considered. Plotted also are the results of the batch reactor (SIR) problem, for the same valueof R ( ρ, . While the results are quite similar, diffusion does affect the shape of the curves obtained, leading to aslightly smaller infection intensity, essentially corresponding to a slightly smaller effective R ( ρ, . Figure 11 is a plotof the asymptotic value r ∞ and of the maximum in infections i max as a function of R ( ρ, for the respective twocases. We find that diffusion acts to moderate somewhat the contagion intensity.We conclude that the two solutions, corresponding to the contagion wave (equations (49)-(50)), now denoted by i D ( t ) .and to the batch (SIR) problem (equations (21)-(22), now denoted by i B ( t ) ), are approximately equal, although notidentical i ( x, t ) → i ( x − V t ) ≡ ˆ i ( ζ ) ≡ i D ( const − xV + t ) ≈ i B ( const − xV + t ) (66)Note that the constant in (66) can absorbed in x/V .The variation of the dimensionless velocity, equation (64) with R ( ρ, is shown in Figure 12. As expected it increaseswith R ( ρ, , varying roughly in a linear fashion at relatively large values of the latter. We conclude that diffusion is14 PREPRINT - A
UGUST
31, 2020Figure 8: Infected and susceptible density fraction profiles at different values of time ( R ( ρ,
0) = 2 . ).fundamental to the spreading of the contagion, the rate of spreading increasing with the square root of the diffusioncoefficient and with the value of R ( ρ, . It also affects, although only mildly, the overall infection profiles, by loweringtheir intensity. Restricting mobility, here represented by D , confirms a most important policy effect for containinginfection. The finding that diffusion also affects somewhat the effective value of R ( ρ, , is a result intuitively expected,but not previously identified. And the fact that the limit of small diffusion is a singular one calls for some caution on theinterpretation of any batch problem results.The wavelike spread of infections predicted by the model has been observed in several pandemics. In this context it isinteresting to compare the spread of the 1347-1350 black death (pneumonic plague) in Europe and that of the 2009pandemic influenza in the US. While there is significant debate regarding the effective R value for these pandemics,it is generally accepted that they are in the range of R ≈ . − . [6, 29]. Further, for both diseases the infectiousperiod is around two weeks, and therefore Λ ≈ / days − . To predict the wave speed, we need an estimate of therespective diffusion coefficients. It is to be expected that in the medieval times, the diffusion of humans, therefore ofinfection, was much smaller. In [6, 30] this was determined by how quickly information was disseminated (in thatcase by word of mouth), resulting into a diffusivity estimate of about miles /day . Using the obtained estimate in15 PREPRINT - A
UGUST
31, 2020Figure 9: Infected density fraction as a function of space and time coordinates ( R ( ρ,
0) = 2 . ).Figure 10: Variation of susceptible ( s ), infected ( i ) and recovered ( r ) fractions (steady state) as a function of non-dimensional time with diffusion included (solid lines) and for the batch reactor (SIR) problem (dashed lines) for thesame value of R ( ρ,
0) = 2 . . The initial infection fraction for the batch system is i B (0) = 0 . our estimate for the speed of the contagion (see Figure 12) we arrive at a speed of . miles/day for the black deathpandemic, which is in the ball-park of the actual observed speed of around mile/day [31].For the 2009 influenza pandemic, we can make a similar back-of-the-envelope calculation. Since children play animportant role in the spread of influenza, we derive this estimate based on their activities. We assume that a typical childwill see 30 other children in school and playgrounds during a typical day. Further, they will travel an average of 10miles to get to these places. In a coarse analogy to an ideal gas this gives a mean-free path of 10 miles and a frequencyof collision of days − , and therefore a diffusivity of around , miles /day . Using these parameters in ourestimate of the speed of the contagion wave, we arrive at a speed of miles/day for the 2009 influenza pandemic.This is also in the ball-park of the observed wave speed of around miles/day for the 2009 influenza pandemic [20].The higher diffusion coefficient for the 2009 influenza pandemic leads to a significantly faster wave speed.16 PREPRINT - A
UGUST
31, 2020Figure 11: Variation of herd immunity ( r ∞ ) and maximum infection fraction ( i max ) as a function of R ( ρ, for thetwo cases of a batch system (dashed lines) and when diffusion is included (solid lines). The initial infection fraction forthe batch system is i B (0) = 0 . .Figure 12: Variation of non-dimensional wave speed in a one-dimensional problem as a function of R ( ρ, .For completeness, we simulated the “collision” of two waves, one moving from the left and the other from the right.Figure 13 shows how the two waves amplify as they interact, then ultimately decay as infection has spread fully and thepopulation reaches conditions of herd immunity at the corresponding value of R . This wave interaction is different, fora number of reasons, from what is observed in other wave problems (e.g. solitary waves, and where for example, thewave velocities increase with the wave amplitude [32]). In the present context, the wave amplitude, e.g. i max cannotexceed the value of 1.We conclude this section by considering effects of heterogeneity in this 1-D rectilinear geometry. Figures 14 and 15show results when an infection wave enters a region with a lower value of R ( ρ, , e.g. one of lower spatial density,from a region of a higher value of R ( ρ, , e.g. one of higher spatial density. For example, such could be the case ofcontagion spreading from an urban to a rural area. In the figure this occurs at x = − . As it enters the low-densityregion the wave decelerates, with the magnitude of infected fraction decreasing. The slowing of the wave is indicated17 PREPRINT - A
UGUST
31, 2020Figure 13: Infected and susceptible density fraction profiles at different values of time for two waves propagatingtowards the center of the domain ( R ( ρ,
0) = 2 . ).by the increase of the slope in the x − t domain. Conversely, at x = 25 , when the wave now enters a region of higher R ( ρ, , it accelerates to a larger velocity, the intensity of the infected fraction correspondingly increasing. The results of the 1-D rectilinear geometry also apply almost identically to radial geometries. Now, equations (49)-(50)become ∂s∂t = Cµ ∂∂µ ( µ ∂s∂µ ) − R ( ρ, r ) si (67) ∂i∂t = Cµ ∂∂µ ( µ ∂i∂µ ) + R ( ρ, r ) si − i (68) ∂r∂t = Cµ ∂∂µ ( µ ∂r∂µ ) + i (69)18 PREPRINT - A
UGUST
31, 2020Figure 14: Infected and susceptible density fraction profiles at different values of time in a 1-D heterogeneous system.For x ∈ ( − , − and x ∈ (25 , , R ( ρ,
0) = 3 , while for x ∈ ( − , , R ( ρ,
0) = 1 . .where we denoted the radial coordinate with µ . We look for the solution of the problem when a region around the origin( µ = 0 ) is initially infected. As in the rectilinear geometry case, we expect that the solution will evolve in terms of atraveling wave, therefore, we consider a transformation to the moving coordinates ξ = µ − V t , and t (cid:48) = t . In thesecoordinates, the equations become ∂s∂t (cid:48) − V ∂s∂ξ = C ∂ s∂ξ + Cξ + V t (cid:48) ∂s∂ξ − R ( ρ, r ) si (70) ∂i∂t (cid:48) − V ∂i∂ξ = C ∂ i∂ξ + Cξ + V t (cid:48) ∂i∂ξ + R ( ρ, r ) si − i (71) ∂r∂t (cid:48) − V ∂r∂ξ = C ∂ r∂ξ + Cξ + V t (cid:48) ∂r∂ξ + i (72)We note that in the limit of large t (cid:48) , equations (70)-(71) transform to the same equations as (52)-(53), with the firstterm on the left and the second term on the right of the equations vanishing. Therefore, the same results hold for19 PREPRINT - A
UGUST
31, 2020Figure 15: Infected density fraction as a function of space and time coordinates in a 1-D heterogeneous system: For x ∈ ( − , − and x ∈ (25 , , R ( ρ,
0) = 3 , while for x ∈ ( − , , R ( ρ,
0) = 1 . .radial geometries as those for the rectilinear problem. The simulations for the 1-D radial geometry, shown in Figure16, confirm the theoretical findings. Interestingly, diffusion in the above reaction systems does not enter in the moreFigure 16: Infected density fraction as a function of radial (space) and time coordinates ( R ( ρ,
0) = 2 . ).familiar form of a similarity variable η = µ √ t , but rather in terms of a wave that propagates with a constant linear speed. Consider, now infection propagation in a general, heterogeneous 2-D system. We are interested in understanding howpropagation occurs and whether or not one can use the asymptotic wave solutions obtained for the 1-D geometries. Inparticular, we are interested in knowing if one can use a wave equation to describe the evolution of the contagion fronts.To explore this question we consider three different geometries, two corresponding to a layered (stratified) system, andone corresponding to a set of four quadrants. 20
PREPRINT - A
UGUST
31, 2020
Layered System
Consider propagation in a stratified layered system, the middle layer having a larger density value(thus, a larger R ( ρ, ), compared to its two adjacent layers. Infection is initiated uniformly on the left boundary (seeFigure 17). As expected, an infection wave emerges in the middle layer, traveling with a velocity v correspondingFigure 17: Infected density fraction ( i ) at t = 40 in a layered medium with an infection wave starting at x = 0 . In theouter layers R ( ρ,
0) = 2 while in the inner layer R ( ρ,
0) = 3 . . The spatial dimension of the domain is × .to the high R value in that layer, and equal to the corresponding asymptotic 1-D wave velocity of 3.15 (see Figure12) for R ( ρ,
0) = 3 . . A slower moving contagion front develops in the adjacent layers, with corresponding velocity v equal to 1.96 which is the asymptotic value of the 1-D wave velocity corresponding to R ( ρ,
0) = 2 . Because ofdiffusion, however, transverse waves also emerge, emanating at the interface of the layers, propagating in the transverse( y ) direction. These waves manifest as a straight line, slanted with respect to the transverse direction, but movingforward with speed v .If we set the origin of the coordinates at the intersection of the interface between the two layers and the initial infectionboundary), the straight line connecting the two fronts in the two layers is described by the following equation F ≡ y ( v − v ) + xa − v at = 0 (73)where the leading edge travels with the middle layer velocity v , the trailing edge with the adjacent layer velocity v ,and the tip of the connecting layer at y = y ( t ) with velocity a , to be determined. Because the waves emanating fromthe interface travel in the adjacent layers at fixed x with transverse velocity v , equation (73) gives v ( v − v ) = v a (74)therefore we find a = v v ( v − v ) (75)The slope of the straight line must then be slope = a ( v − v ) = v v (76)The simulations of (17) show that this is indeed the case.More generally, if we define a front by F ( x , t ) = 0 where x is the space vector, its evolution satisfies F t ( x , t ) + v · ∇ F = 0 (77) F t ( x , t ) + v n |∇ F | = 0 (78)where subscript t denotes time derivative, and v n is its component in the direction of ∇ F . For rectilinear or radialgeometries, as well as in each of the two layers discussed, (78) reduces to a linear wave propagating at a constant21 PREPRINT - A
UGUST
31, 2020velocity. However, for the straight line connecting the two waves, resulting from a sequence of waves emanating fromthe layer interface, the relevant equation is (73), which leads to a normal velocity equal to v n = v (cid:113) (1 + v v ) (79)The interdependence between the three velocities (leading front, v , normal to the surface, v , and trailing front, v ) is tobe noted. Four heterogeneous quadrants
A different illustration of these effects is shown in Figure 18, where infectioninitiates in the upper right corner of a rectangular domain with four different quadrants, each taken with a differentdensity (hence different values of R ( ρ, ). We have selected the two diagonally opposite quadrants (e.g. NE and SW)to have relatively larger values of R ( ρ,
0) = 3 . , hence an asymptotic wave velocity of 3.15, while the other two (SE,NW) to be of relatively smaller values of R ( ρ,
0) = 2 , hence an asymptotic wave velocity of about 2. The infectionwaves follow the expected pattern. Infection grows first radially in the NE quadrant with the asymptotic speed of 3.15(upper right panel in Figure 18). When the SE and NW boundaries are encountered (lower right panel in Figure 18),the waves slow down and start spreading in a radial manner as they enter the two quadrants at the slower speed of 2.Subsequently, the infection wave enters the high density SW region, in which it starts spreading radially, now movingwith the higher velocity of 3.15 (lower left panel in Figure 18). Upon touching the boundaries with the NW and SEregions, however, diffusion causes infection waves to start emanating from the higher to the lower density regions,resulting into a straight line segment, similar to that for the layered system, connecting the two different waves in thetwo different regions. For the same reasoning as in the layered system, this straight line segment has a slope equals tothe ratio of the two respective velocities, namely (76). The simulations confirm these results.Figure 18: Contour plot of the infected density fraction a region with four quadrants with different values of R ( ρ, atdifferent values of time. R ( ρ,
0) = 3 . for NE and SW, and R ( ρ,
0) = 2 for NW and SE, respectively. The spatialdimension of the region is × .We end this section by noting that one can further explore the use of a wave equation, or perhaps the eikonal equation,for the description of the spreading of contagion. This would be the goal of a separate study.22 PREPRINT - A
UGUST
31, 2020
We used an analogy with chemical reaction engineering processes to model the growth and spreading of human-transmitted infections, such as COVID-19 and influenza. The basis of the model is the assumption of three distinctpopulations, as in the celebrated SIR model, which are mapped into equivalent chemical species. An important firstresult from this formulation is that the relevant quantities are population densities (specifically, areal spatial densities)and not total populations, as is in the SIR model. These satisfy partial differential equations, with spatial mobilitiesincluded through diffusion and advection. Spatial density effects are then incorporated into the kinetic constant ofthe infection growth, hence into the key parameter R , which becomes an explicitly function of spatial densities, andfound to decrease with the extent of the contagion. Incorporating a spatial density dependence in contagion kinetics isnecessary, and consistent with health guidelines and droplet dynamics.In the absence of diffusion, and assuming spatially uniform profiles, the results revert to a modified version of the SIRmodel, in a geometry recognized as a “batch reactor”. We find that if were to identify an effective SIR-like R , it wouldbe smaller than when spatial effects are ignored, which is the standard SIR case. The infection curves are found to belargely independent of initial conditions, the effect of which is simply to delay the onset of the infection. Analogousresults hold when a ban on “imported infection” is applied, which also only affects to delay the onset of the epidemic,assuming that R is not modified. Small fluctuations in the initial conditions have a minimal effect on the ensemblebehavior. However, this is not the case when the fluctuations are significant in space or in time, in which case they resultinto non-trivial ensemble averages. Using the density-dependent R we can readily model the effect of the equivalentof a commute between “home” and “work”, or “suburban” and “urban”, to obtain an effective R , which is found to beequal to the arithmetic mean weighted by the exposure time.We subsequently considered the effect of diffusion in 1-D geometries. We showed the emergence of infection waves,whether in rectilinear or radial geometries, which asymptotically travel with constant velocity, the dimensional valueof which scales with the square root of diffusivity, and increases with R . The behavior of the infection waves ata fixed time is very similar to that for an SIR system, although the shape and the effective value of R are slightlydifferent than in the absence of diffusion, leading to lower intensities, than in the absence of diffusion. This result isimportant, as it shows that the SIR limit (absence of diffusion) is a singular limit. We then examined how distributeddensities in different geometries affect the propagation of infection waves. We find that for all practical purposes, thecontagion propagation can be modeled by equivalent waves, the speed of which only depends on the value of R in thatregion. The overall conclusion is that spatial effects via density or diffusion considerations are important and must beconsidered in the description of contagion and its epidemics.23 PREPRINT - A
UGUST
31, 2020
References [1] William Ogilvy Kermack and Anderson G McKendrick. A contribution to the mathematical theory of epidemics.
Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character ,115(772):700–721, 1927.[2] Fabrizio Croccolo and H Eduardo Roman. Spreading of infections on random graphs: A percolation-type modelfor covid-19.
Chaos, Solitons & Fractals , page 110077, 2020.[3] Tiberiu Harko, Francisco SN Lobo, and MK Mak. Exact analytical solutions of the susceptible-infected-recovered(sir) epidemic model and of the sir model with equal death and birth rates.
Applied Mathematics and Computation ,236:184–194, 2014.[4] Roy M Anderson and Robert M May. Population biology of infectious diseases: Part i.
Nature , 280(5721):361–367,1979.[5] Keith Burghardt and Kristina Lerman. Unequal impact and spatial aggregation distort covid-19 growth rates. arXiv preprint arXiv:2004.12994 , 2020.[6] Julia V Noble. Geographic and temporal development of plagues.
Nature , 250(5469):726–729, 1974.[7] Sheldon Ross.
A First Course in Probability 8th Edition . Pearson, 2009.[8] Felix Sharipov and José L Strapasson. Ab initio simulation of transport phenomena in rarefied gases.
PhysicalReview E , 86(3):031130, 2012.[9] H Scott Fogler.
Essentials of Chemical Reaction Engineering: Essenti Chemica Reactio Engi . Pearson Education,2010.[10] R Byron Bird. We stewart, and en lightfoot.
Transport phenomena , 11:5, 1960.[11] Hydrodynamics Lamb and Sixth Ed Hydrodynamics. Cambridge university press, 1932. et seq , page 254, 1966.[12] Ronald R Coifman and Stéphane Lafon. Diffusion maps.
Applied and computational harmonic analysis ,21(1):5–30, 2006.[13] Franca Hoffmann, Bamdad Hosseini, Assad A Oberai, and Andrew M Stuart. Spectral analysis of weightedlaplacians arising in data clustering. arXiv preprint arXiv:1909.06389 , 2019.[14] Benoit B Mandelbrot. The fractal geometry of.
Nature , pages 394–397, 1982.[15] AS Chaves. A fractional diffusion equation to describe lévy flights.
Physics Letters A , 239(1-2):13–16, 1998.[16] Diego del Castillo-Negrete, BA Carreras, and VE Lynch. Front dynamics in reaction-diffusion systems with levyflights: a fractional diffusion approach.
Physical Review Letters , 91(1):018302, 2003.[17] Péter Érdi and János Tóth.
Mathematical models of chemical reactions: theory and applications of deterministicand stochastic models . Manchester University Press, 1989.[18] Qun Liu and Daqing Jiang. Dynamical behavior of a stochastic multigroup sir epidemic model.
Physica A:Statistical Mechanics and its Applications , 526:120975, 2019.[19] JA Wesselingh, Rajamani Krishna, et al.
Mass transfer in multicomponent mixtures . Delft University Press Delft,2000.[20] Julia R Gog, Sébastien Ballesteros, Cécile Viboud, Lone Simonsen, Ottar N Bjornstad, Jeffrey Shaman, Dennis LChao, Farid Khan, and Bryan T Grenfell. Spatial transmission of 2009 pandemic influenza in the us.
PLoS ComputBiol , 10(6):e1003635, 2014.[21] Alex Viguerie, Guillermo Lorenzo, Ferdinando Auricchio, Davide Baroli, Thomas JR Hughes, Alessia Patton,Alessandro Reali, Thomas E Yankeelov, and Alessandro Veneziani. Simulating the spread of covid-19 viaspatially-resolved susceptible-exposed-infected-recovered-deceased (seird) model with heterogeneous diffusion. arXiv preprint arXiv:2005.05320 , 2020.[22] Joseph O Hirschfelder, Charles F Curtiss, Robert Byron Bird, and Maria Goeppert Mayer.
Molecular theory ofgases and liquids , volume 165. Wiley New York, 1964.[23] William F Wells. On air-borne infection: Study ii. droplets and droplet nuclei.
American journal of Epidemiology ,20(3):611–618, 1934.[24] Lydia Bourouiba, Eline Dehandschoewercker, and John WM Bush. Violent expiratory events: on coughing andsneezing.
Journal of Fluid Mechanics , 745:537–563, 2014.[25] X Xie, Y Li, ATY Chwang, PL Ho, and WH Seto. How far droplets can move in indoor environments–revisitingthe wells evaporation–falling curve.
Indoor air , 17(3):211–225, 2007.24
PREPRINT - A
UGUST
31, 2020[26] Vasilis Marmarelis. Predictive modeling of covid-19 data in the us: Adaptive phase-space approach.
IEEE OpenJournal of Engineering in Medicine and Biology , 2020.[27] Athanasios S Fokas, Nikolaos Dikaios, and George A Kastis. Covid-19: Predictive mathematical models for thenumber of deaths in south korea, italy, spain, france, uk, germany, and usa. medRxiv , 2020.[28] Athanasios S Fokas, Jesús Cuevas-Maraver, and Panayotis G Kevrekidis. Two alternative scenarios for easingcovid-19 lockdown measures: one reasonable and one catastrophic. medRxiv , 2020.[29] Roya Nikbakht, Mohammad Reza Baneshi, Abbas Bahrampour, and Abolfazl Hosseinnataj. Comparison ofmethods to estimate basic reproduction number (r0) of influenza, using canada 2009 and 2017-18 a (h1n1) data.
Journal of research in medical sciences: the official journal of Isfahan University of Medical Sciences , 24, 2019.[30] Nicolas Rashevsky.
Looking at history through mathematics . Cambridge, Mass: MIT Press, 1968.[31] William L Langer. The black death.
Scientific American , 210(2):114–121, 1964.[32] Philip G Drazin and Robin S Johnson.