key: cord-0524652-lcnzmc3v
authors: Vyasarayani, C. P.; Chatterjee, Anindya
title: Complete dimensional collapse in the continuum limit of a delayed SEIQR network model with separable distributed infectivity
date: 2020-04-26
journal: nan
DOI: nan
sha: dcd258e63a7c3af8e305c4154433995bd82a6ba1
doc_id: 524652
cord_uid: lcnzmc3v

We take up a recently proposed compartmental SEIQR model with delays, ignore loss of immunity in the context of a fast pandemic, extend the model to a network structured on infectivity, and consider the continuum limit of the same with a simple separable interaction model for the infectivities $beta$. Numerical simulations show that the evolving dynamics of the network is effectively captured by a single scalar function of time, regardless of the distribution of $beta$ in the population. The continuum limit of the network model allows a simple derivation of the simpler model, which is a single scalar delay differential equation (DDE), wherein the variation in $beta$ appears through an integral closely related to the moment generating function of $u=sqrt{beta}$. If the first few moments of $u$ exist, the governing DDE can be expanded in a series that shows a direct correspondence with the original compartmental DDE with a single $beta$. Even otherwise, the new scalar DDE can be solved using either numerical integration over $u$ at each time step, or with the analytical integral if available in some useful form. Our work provides a new academic example of complete dimensional collapse, ties up an underlying continuum model for a pandemic with a simpler-seeming compartmental model, and will hopefully lead to new analysis of continuum models for epidemics.

The global pandemic of COVID-19 has prompted several studies of epidemic models from a dynamic systems point of view. Pandemics can be studied using simple mean-field models or compartmental models, where the entire population is divided into susceptible (S), exposed (E), infectious (I), quarantined (Q), and recovered (R) groups. More groups like hospitalized (H), vaccinated (V), etc., can be added, or groups can be removed, depending on modeling goals. These models are primarily developed based on the original SIR model due to Kermack and McKendrick [18] . Many models include other complexities like prior immunity, temporary immunity transferred at birth, vaccination history, a carrier population that never recovers [14] , reinfection due to loss of immunity after recovery [6] , exposed but asymptomatic populations [13] , a quarantined population [9] , and the influence of vital dynamics [12] .

The fidelity of such models can be improved by developing structured network models [29, 20, 11, 26, 7, 17] , where each node in the network is a compartmental model with different parameters. In other words, the entire population is divided into N subgroups or structures. Each such subgroup or structure can have its own S, E, I, Q, and R subpopulations. These subgroups or structures can be based on age, lifestyle, demography, geography, or other aspects. The network topology [17] can be uniform, timevarying, or random. Recently, for example, the initial spread of COVID-19 in India was explored using an age-structured SIR network model with simulated social lockdown conditions [25] .

In contrast to such network models, an alternative high-dimensional approach uses reaction-diffusion [2, 23, 16, 24] type partial differential equations (PDEs) to capture the spatial and temporal evolution of an epidemic. Such models with spatial dimensions allow solutions with concentrated pockets of infection which emerge and spread. Sometimes, PDEs with a single spatial variable are converted into integro-differential equations [19, 23] .

Time delays [8] are often present in the progression of diseases due to latency and incubation times. Several researchers have studied the spread of infections like Zika, HIV, Hepatitis, and Influenza using delay differential equations [28, 5, 21, 1, 10, 22] . Recently, an SEIQR [28] model with time delays was proposed to study the progression of a generic epidemic. In recent work of our own [27] , we have studied the time delayed SEIQR model of [28] after neglecting loss of immunity over time, which is appropriate for a fast pandemic.

In the lumped approach of [28, 27] , the infectivity or transmission rate of the disease is modeled using a single positive parameter β =βm.

In Eq. 1,β is a characteristic of the pathogen; and m, the density of contacts, is a characteristic of the behavior of individuals.

In this paper, we address the situation where the parameter β is a distributed quantity over the entire population. After all, some people have greater exposure to infection than the population average (e.g., police personnel, people providing other essential services, hospital staff, as well as people less willing to cooperate with government-recommended social distancing measures). Others have less exposure than average (e.g., older people, people who cooperate more with recommended social distancing measures). Some have strong immune systems, and others have weak immune systems, which might be reflected in differing values of the underlying parameterβ in Eq. (1). All these people interact to various degrees in modern society, where travel and mixing are common. The details of this variation, if incorporated, will yield a richer model that can hopefully make more realistic predictions. We clarify that we do not stratify the population by age, occupation, health status, or social behavior. In our model, all these factors contribute to a final effective β for each person. We stratify the population according to β alone, because that is what affects the dynamics of infection.

We ask the following question: for a general distribution of β in the population, and for a reasonable model for how different sections with different β's interact during the pandemic, how does the continuum solution evolve from infinitesimal initial infection all the way to final saturation?

In the rest of this paper, we show the following. Under a simple separable interaction model for different β's, and upon neglecting loss of immunity for those who have recovered from the disease, the continuum model reduces to two nonlinear integrodifferential equations with delays, with the continuously variable β appearing as a parameter. However, due to a remarkable dimensional collapse, the dynamics is exactly described by a single scalar nonlinear delay differential equation without integrals, where the parameter β appears only through the first derivative of its moment generating function. If some required moments of the distribution are finite, then a local expansion can be easily computed for small levels of overall infection, and the correpondence with the lumped model of [28, 27] is transparent and close. If the moments are unbounded, alternative, slightly more complex, single DDE approximations can still be developed, at least in principle. In this way, both the justification for the lumped model, as well as corrections needed for variable β, are exactly demonstrated. From a practical viewpoint, our results clarify the role of, e.g., higher versus lower variability of β in the population in the spread, saturation, and containment of the disease. From an academic viewpoint, our work presents a new and satisfying example of extreme model order 2 reduction.

In [28, 27] , a single fixed parameter β is used to model the whole population. Our presentation begins with a brief statement of this model. Figure 1 shows the schematic of the single-β SEIQR model and is adapted from [28] . The governing equations, included for completeness, are: Figure 1 : The SEIQR model with delays for each node. Here, S(t) represents healthy individuals. The infection rate constant is β. Asymptomatic and infected individuals E(t) remain non-infectious for σ units of time. Later, they become infectious and are represented by I(t), but show no symptoms for another τ units of time. When symptoms appear, these infected people are isolated or quarantined with probability p for a time κ, and are represented by Q(t). A few asymptomatic but infectious individuals may recover on their own, at a rate γ. The cured population R(t) after quarantine may lose immunity at a small rate α, but we use α = 0 over the time scale of a fast-spreading pandemic.

In Eq. (2), the right hand side has two terms: the rate of change of the susceptible population S(t) depends on an infection rate βS(t)I(t) and a rate of reinitroduction of susceptible people through loss of immunity modeled using αR(t). When β varies within the population, we use a network model. Figure 2 illustrates the idea. Here the total population is partitioned into N sub groups P 1 ,P 2 ,...,P N ; and each group is modeled with different β k ranging from β 1 to β N . The primary consequence is that in the network's equivalent of Eq. 2 for group k, the I(t) must be replaced by something that represents interaction effects from all groups. We will model this using a quantity Λ(t) that applies to the entire network, as explained below. Figure 2 : A schematic representation of N interacting population groups with different infection spread rates among each group. Every connection between two groups is bidirectional and symmetric, and every group is connected to all other groups (a dense network).

The main idea is simple. A susceptible person in group k, for any 1 ≤ k ≤ N , can get infected through interaction with a person in any group r, with 1 ≤ r ≤ N . Thus, the total infection rate of the single-β model, i.e., βS(t)I(t), is to be replaced by the sum of contributions from all groups. The contribution from one single group r is taken to be

for some suitable interaction coefficient F (β k , β r ) that has to be specified as part of the mathematical model. The function F (β k , β r ) is understood to be symmetric [4] , i.e.,

It must monotonically increase with respect to each of its arguments, i.e.,

Clearly, we require F (x, 0) = 0 regardless of x, because a person who meets nobody at all will not get infected. Further, it seems reasonable to assume that if a person who mixes very little meets another person who mixes very little, then the probability of infection spreading is much lower than if the second person mixes widely. Mathematically, we assume that

All the above conditions are met by the separable choice

Finally, in the special case where the entire population moderates its behavior to produce a single effective β, the net effect should ideally reduce to Eq. (2) , and this is something that will turn out to be true for Eq. (8).

The above choice of F (β k , β r ) leads us to define

and the dynamics of the network is governed bẏ

In this model, except for β k , all other parameters namely σ, τ , κ, p, γ, and α are the same for each node in the network. This assumption is motivated by the idea that, in a well-functioning society, the degree of mixing practiced by individuals varies a lot more than the time taken for an infected person to be detected and quarantined; and that the rate of self-recovery γ is a biological quantity indpendent of an individual's social behavior. Note that we assume the loss of immunity rate, α, to be zero. We observe from Eqs. (10)-(14) that if α = 0 the states E k (t), Q k (t), and R k (t) become slave variables. Only S k (t) and I k (t) need to be solved for. Further, we can always scale time by setting σ = 1, making τ , β k , and γ effectively dimensionless. Equations (10)- (14) , after droppingĖ,Q andṘ and defining

We have integrated Eqs. (15) and (16) using Matlab's built-in solver dde23 for many different initial functions or history functions, with error tolerances set to be 10 −7 or better. For ease of presenting results, we define the arraysβ

The history functions for numerical integration are selected in terms of nonnegative functionsÛ (β) and

In the above,Û (β) was chosen to satisfy ||Ŝ(β, −ν)|| 1 = 1. It is clear that ||Î(β, −ν)|| 1 = 0 regardless of V (β). Such initial conditions, with small initial infection, are of interest because we wish to study whether the infection remains small or grows significantly within the population. It is clear thatÛ (β) plays the role of a probability density function of people's β-values in the population. IfÛ (β) is large for small β and decays rapidly for large β, then most people are in a regime of low infection rates. Conversely, ifÛ (β) decays slowly with increasing β, then a signifcant proportion of the population is in a regime of high infection rates. The dynamic consequences of such different distributions will be studied using initial numerical simulations in this section, before proceeding to analytical treatment. The net effect of the pandemic, or net damage, is taken to bê

For total percentage affected, we use

For our initial numerical case study, we use N = 500 and select β values uniformly spaced between β 1 = 0 and β N = 6, i.e., β values are distributed within a strictly finite range. Figure 3 (a) showsŜ(β, t) at three different time instants at t = −ν (initial time), t = 40, and t = 100. Also shown is a fit of the form

at these three instants of time. The initial conditions werê

which matches the Gamma distribution, with a = 5 and b = 0.4 (technically, the Gamma distribution has infinite support, but for the parameters chosen it has decayed to tiny values for β = 6). Also, the other initial functionV (β) was taken to be identical toÛ (β) in the simulations of figure 3. We emphasize that in Fig 3(a) , f (−ν) = 0 by definition and there are only two fitted numbers: f (40) and f (100). The match is essentially perfect, and shows that the evolution of all the S k (t) together have a one-dimensional behavior. The solution at any time t is fitted, essentially perfectly, by a function of the formŜ

where "•" denotes elementwise multiplication of two arrays. Figure 3 √ β is essentially exact. In other words, the variation over β is one-dimensional, in terms of a scalar f (t). We refer to this great reduction in dimensionality, where the variation with respect to the continuous variable β is accounted for by a single scalar, as complete dimensional collapse. Numerics indicate the collapse is exact.

It remains only to extract the governing equation for the scalar f (t), and we will do this using a continuum formulation.

We consider a continuum limit of Eqs. (15) and (16) , as N → ∞. We replace the summation in Eq. (9) with integrals and assume S(β, t) and I(β, t) to be functions of β and t. In the continuum limit, Eqs. (15) and (16) becomeṠ

In the above integrals, if β takes only finite values (in the numerical examples above, β was between 0 and 6), then the upper limits of the integrals can be replaced with finite values; we write ∞ as a formal upper limit. Note that Eqs. (20) and (21) represent infinitely many coupled delay differential equations, parameterized by the distributed infectivity parameter β.

The fraction of the population that is susceptible is now understood to be ∞ 0 S(β, t) dβ, and the fraction of the presently infectious population is understood to be

At the start of the pandemic, we have initial conditions that satisfy Here, finitely supported polynomials of the formÛ (β) = r n · 1 +β n , 0 ≤β ≤ 6 are considered, with r n chosen to make the sum (or 1norm) ofÛ (β) equal to 1. The distribution of initial values forÎ(β) usesV (β) that is random, initially uniformly distributed between 1 and 2, and then normalized to unit 1-norm. Also shown is the fitŜ(β, t) = S(β, −ν)e −f (t) √ β at different instants of time. These results demonstrate that provided the intial infected population distribution is small, the subsequent evolution obeys a simple one-dimensional description with a scalar variable f (t). 9 has greater social mixing, i.e., S(β, 0) that decays slowly with increasing β, may see an outbreak of the infection. A point to note during numerical solution of Eqs. (20) and (21) is that S(β, t) and I(β, t) are nonnegative. It is in principle possible for I(β, t) to start positive and become exactly zero at some instant (based on hypothetical initial functions used in the delayed variables), and it corresponds to all infected people becoming quarantined before fresh people are infected. However, such a situation corresponds to the pandemic being quenched by eliminating infection, and is not of interest when either infection cannot be eliminated, or when even elimination of infection leaves the system unstable (i.e., introduction of infinitesimal infection leads to an outbreak). As seen in the numerical solutions of Figs 3 and 4 , there are clearly also other solutions of interest where the pandemic progresses, infects a percentage of the population, and reaches an eventual equilibrium with the infection not progressing further. For such solutions, the nonnegativity constraint of S(β, t) and I(β, t) turns out to be satisfied. 33) is taken as f (t) = ω × (1 + t ν ), with ω is adjusted to match the results from (Eqs (15) and (16)). Other numerical parameters are mentioned in text boxes within subplots.

Motivated by the foregoing observations from numerical simulations, we now look for solutions of the form

where φ(β) can be determined from initial conditions. Substituting into Eq. (20) yields

The integral on the right hand side is purely a function of time t: let us call it g(t). Then we havė

Also, Eq. (21) becomesİ

We multiply both sides above by √ β and integrate with respect to β. Defining

and assuming the integral and the derivative on the left hand side can be interchanged, Eq (25) becomeṡ

Any solution of Eqs. (24) and (27), with H(t) defined by Eq. (26) , corresponds to an exact solution of Eqs (20) and (21) . Differentiating Eq. (24) with respect to time and substituting forġ(t) from Eq. (27), we obtain:f

Separately, multiplying both sides Eq. (26) byḟ (t), we obtain

where we notice that the right hand side is integrable with respect to time, giving

where C is a time-independent constant (it does not depend on β either, because β here is a dummy variable of integration). Let the integral on the right hand side of Eq (30) be called G(f (t)), i.e.,

where we notice that the time t within this definition is merely a parameter that determines f , and so we can also simply write

if it suits our purpose. Equation (28) becomes, upon one integration,

where C 0 is a constant of integration and can be evaluated by setting f (t) = 0 whenḟ (t) = 0 (at the time of infinitesimal intiation of infection), and is given by

For many initial distribution functions φ(β), G(f (t)) and C 0 can be evaluated in closed form. For example, for the uniform distribution

where

In this way, at least in principle (even if the integrals cannot be evaluated in closed form), Eq. (33) is a first order scalar DDE that governs the dynamics of the continuum limit of our network. The infinitedimensional variability with respect to the distributed parameter β has collapsed into a single dimension. This observation of complete dimensional collapse is one of the main contributions of the paper. Our result is independent of the distribution φ(β). Figure 5 compares the continuum solution −f (t) obtained from Eq. (33), and the finite network solution (Eqs (15) and (16)

β for different parameter values (see figure caption for details). The results match closely.

We now proceed to interpret G(f ) in terms of the moment generating function of an underlying randomly distributed quantity,

by defining an intermediate quantity (using an overbar to formally distinguish the two functions)

For example, if φ(β) = β, then we meanφ(u) = u 2 . Equation (31), written as Eq. (32), becomes

Now, if we interpret φ(β) as the probability density function of the random variable β in the population, and think of ψ(u) as the probability density function of the transformed variable u = √ β in the same population, then it must be true that

which gives

where we see that, except for a sign change in f , G(f ) is the first derivative with respect to f of the moment generating function [15] of the distributed quantity u = √ β. If we expand in a series for small f , we obtain

where the coefficients m k are moments of u, given by

assuming these moments exist. In other words, m 1 is the population mean of √ β, m 2 is the population mean of β, m 3 is the population mean of β 3/2 , and so on. Returning to Eq. (33) and expanding for small f , and retaining up to cubic terms (i.e., retaining up to second order nonlinear correction terms), we havė

In Eq. (42), if we are interested in an outbreak that starts from an infinitesimal infection, such that f = 0 whenḟ = 0, then the constant

consistent with Eq. (34), leavinġ

Linearising Eq. (43) for small f , we obtaiṅ

Substituting f (t) = e st , we find s = m 2 e −s − m 2p e −νs − γ

This characteristic equation matches one studied in [28, 27] , and the stability condition is known to be (recall Eq. (22))

In fact, since m 2 is the second moment of √ β, i.e., the population average of β, this validates the use of an average value in the lumped model (Eqs. (2)-(6)), studied in [28, 27] . Since our continuum model allows a general distribution for β, we can consider the specific case of a Dirac delta function (i.e., β is not random)

Equation (40) yields

Equation (33) becomes (inserting C 0 to match initial conditions of f = 0 whenḟ = 0)

Setting f (t) = √ β 0 P (t), we obtaiṅ

which exactly matches Eq. (15) of [27] . Further, even for variable β, the correpondence between Eq. (33) and Eq. (49) actually holds up to second order. If the expansion for small f in Eq. (43) is truncated at second order, i.e., we retain only the second and third moments of u, then it is easy to show that a simple scaling of f makes that equation match the second order expansion of Eq. (49). The implication is that, for relatively small outbreaks, the dynamics under distributed β (i) depends on the expected values of β and β 3/2 , and (ii) is the same, up to a linear scaling, as the dynamics with a single fixed β as studied in [27] . For larger outbreaks, higher moments of the distribution of β begin to play a role, and the match with Eq. (49) deteriorates. , with "numerical" and "analytical" implying the same as in (a). The initial function used was f (t) = 1×10 −9 +1×10 −9 1 + t ν , t ≤ 0.

Finally, we close with an example where the required moments in Eq. (43) do not all exist. Let the probability density function of β be

for which, in Eq. (43), m 2 = 2a and m 3 = ∞. This distribution has a long tail. The integral G(f ) for f > 0 (recall Eq. (40)) can be found (using Maple) in terms of special functions, and the first few terms of an asymptotic series for small f > 0 is as shown below: Fig 6. The match is essentially perfect, which indicates that finiteness of moments of u, and hence the moments of β, is not needed for the underlying DDE to have useful solutions.

In this paper we have considered the extension of a compartmental model for an epidemic or pandemic to a network model, structured on the infectivity parameter β rather than background factors like age or occupation. We have considered a simple interaction model between susceptible individuals with different infectivity parameters β k and β r , and examined the case where the effective interaction parameter between these groups is monotonic in each parameter, symmetric, and separable. Under reasonable qualitative requirements, we settled on an interaction term of the form √ β k β r . With this interaction, the continuum limit of the network model yielded a pair of delayed integrodifferential equations parameterized by β. However, numerical solutions also suggested a remarkable dimensional collapse in the underlying dynamics, such that the entire variation due to the continuously distributed β could be captured using a single function of time. The dynamics in terms of this function, f (t), turned out to be given by a single first order nonlinear DDE. We have referred to this dramatic simplification as complete dimensional collapse.

We are aware of similar dimensional collapse for the Fokker-Planck equation for a plasma in an ion trap (a parametrically forced system without delays): see [3] and references therein. However, we are not aware of such collapse being previously noted for any models within mathematical epidemiology.

The DDE governing the new variable f (t) has strong connections with the lumped-model DDE where β is assumed to be a single constant for the entire population. When the first few moments of the underlying parameter u = √ β in the population are finite, a direct correspondence up to second order can be established with the single-β equation. Even if the required moments are not finite, the DDE remains well posed, although its small-f behavior may be more complicated.

There are interesting social implications that emerge from our model. For example, in the finitemoments case, even if β is distributed such that a small proportion of the population has rather high β values, there is no immediate catastrophic growth because that subpopulation, being small, tends to interact mostly with other parts of the population whose β values are smaller. Even in a population where some people engage in behavior corresponding to high β (i.e., not obeying social distancing), low-β behavior from other members of the population can help to both limit the pandemic and provide differential benefits based on adopted risk levels. Satisfyingly, the criterion for unstable growth from small initial infection is governed solely by the population average of β, although the final saturation value of the infection depends on higher moments. In particular, for a small outbreak due to weak instability under infinitesimal initiation, the saturation value is essentially determined by the population averages of β and β 3/2 , or of u 2 and u 3 respectively.

Future work with such continuum models may incorporate additional features for greater realism, including non-separable interaction terms, behavior modification by the population as the outbreak progresses, new approximate solutions for various special cases including cases with unbounded moments of u, and direct comparisons of model predictions with national or international data.

A delay differential model for pandemic influenza with antiviral treatment

Spread of infectious diseases in a hyperbolic reaction-diffusion susceptible-infected-removed model

Unifying averaged dynamics of the fokker-planck equation for paul traps

An introduction to compartmental modeling for the budding infectious disease modeler

Numerical modelling in biosciences using delay differential equations

Mathematical models in population biology and epidemiology

Epigrass: a tool to study disease spread in complex networks

Time delay in epidemic models

An SEIQR model for childhood diseases

Dynamics of a delay differential equation model of Hepatitis B virus infection

Efficiency of prompt quarantine measures on a susceptible-infectedremoved model in networks

Three basic epidemiological models

The mathematics of infectious diseases

The basic epidemiology models: models, expressions for R0, parameter estimation, and applications

Introduction to mathematical statistics

Differential equations and mathematical biology

Networks and epidemic models

A contribution to the mathematical theory of epidemics

Spreading disease: integro-differential equations old and new

Six susceptible-infected-susceptible models on scale-free networks

Mathematical analysis of delay differential equation models of HIV-1 infection

A fractional-order model for Zika virus infection with multiple delays

Spatial-temporal dynamics in nonlocal epidemiological models

Modelling SIR-type epidemics by ODEs, PDEs, difference equations and cellular automata-A comparative study

Age-structured impact of social distancing on the COVID-19 epidemic in India

Rapid decay in the relative efficiency of quarantine to halt epidemics in networks

New approximations, and policy implications, from a delayed dynamic model of a fast pandemic

Consequences of delays and imperfect implementation of isolation in epidemic control

Epidemic model with isolation in multilayer networks

We thank Sankalp Tiwari for commenting on this work.