Complete dimensional collapse in the continuum limit of a delayed SEIQR network model with separable distributed infectivity.

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 β . 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 β 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 β appears through an integral closely related to the moment generating function of u = β . 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 β . 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. © Springer Nature B.V. 2020.

Introduction

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 [1]. Many models include other complexities like prior immunity, temporary immunity transferred at birth, vaccination history, a carrier population that never recovers [2], reinfection due to loss of immunity after recovery [3], exposed but asymptomatic populations [4], a quarantined population [5], and the influence of vital dynamics  [6].

The fidelity of such models can be improved by developing structured network models [7-12], 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 [12] can be uniform, time-varying, 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 [13].

In contrast to such network models, an alternative high-dimensional approach uses reaction–diffusion [14-17] 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 [15, 18].

Time delays [19] 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 [20-25]. Recently, an SEIQR [20] model with time delays was proposed to study the progression of a generic epidemic. In recent work of our own [26], we have studied the time delayed SEIQR model of [20] after neglecting loss of immunity over time, which is appropriate for a fast pandemic.

In the lumped approach of [20, 26], the infectivity or transmission rate of the disease is modeled using a single positive parameterIn 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 correspondence with the lumped model of [20, 26] 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 reduction. We mention that a preprint of this article has been uploaded at [27] (https://arxiv.org/abs/2004.12405).

Network model

In [20, 26], 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 [20]. The governing equations, included for completeness, are:

Fig. 1

The SEIQR compartmental model with delays

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 over the time scale of a fast-spreading pandemic, like COVID-19.

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 and a rate of reintroduction of susceptible people through loss of immunity modeled using . 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 ,,...,; and each group is modeled with different ranging from to . 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 that applies to the entire network, as explained below.

Fig. 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 , can get infected through interaction with a person in any group r, with . Thus, the total infection rate of the single- model, i.e., , is to be replaced by the sum of contributions from all groups. The contribution from one single group r is taken to befor some suitable interaction coefficient that has to be specified as part of the mathematical model. The function is understood to be symmetric [28], i.e.,It must monotonically increase with respect to each of its arguments, i.e.,Clearly, we requireregardless 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 thatAll the above conditions are met by the separable choice There could be other formulas for satisfying the conditions mentioned above, but the separable nature of Eq. (8) makes it simple as well as appealing.

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 leads us to defineand the dynamics of the network is governed byIn this model, except for , 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 independent 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 , the states , , and become slave variables. Only and need to be solved for. Further, we can always scale time by setting , making , , and effectively dimensionless. Equations (10)–(14), after dropping , and and definingbecome:

Initial numerical observations

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 or better. For ease of presenting results, we defineIt should be noted that . The history functions for numerical integration are selected in terms of nonnegative functions and asIn the above, was chosen to satisfy . It is clear that regardless of . 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 values of and decays rapidly for large values of , then most people are in a regime of low infection rates. Conversely, if decays slowly with increasing values of , then a significant 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 beFor total percentage affected, we useFor our initial numerical case study, we use and select values uniformly spaced between and , i.e., values are distributed within a strictly finite range. Figure 3a shows at three different time instants at (initial time), , and . Also shown is a fit of the format these three instants of time. Where “” denotes elementwise multiplication of two arrays. The initial conditions werewhich matches the Gamma distribution, with and (technically, the Gamma distribution has infinite support, but for the parameters chosen it has decayed to tiny values for ). Also, the other initial function was taken to be identical to in the simulations of Fig. 3.

Fig. 3

Simulations with different initial values . Here, a Gamma distribution is used for , along with . Also shown is the fit of the form of Eq. (19) at different instants of time. We note that as the values in the distribution get smaller on average, decreases

We emphasize that in Fig 3a, 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 together have a one-dimensional behavior. The solution at any time t is fitted, essentially perfectly, by a function of the form shown in Eq. (19).

Figure 3b shows and the total infected population D(t) during the progression of the pandemic. Parameters used for generating the results are shown in the text boxes inside the figure.

Other simulations with different a and b in the underlying Gamma distribution are shown in subplot pairs Fig. 3c, d and Fig. 3e, f. The quality of the match remains excellent, as is shown graphically for different intermediate time instants in Fig. 3c and e, although the net percentage of people infected changes with a and b (more people are infected if typical values in the distribution are higher).

In Fig 4, we show results for initial conditions chosen arbitrarily on the finite interval [0, 6], without an underlying asymptotic approach to zero for large (see figure caption). The choices for are now polynomials on the interval [0, 6]. From Fig. 4, we observe that for several different and and for generic small-infection history functions [Eqs. (17) and (18)], a fit of the form given by Eq. (19) 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.

Fig. 4

Simulations with different initial values . Here, finitely supported polynomials of the form are considered, with chosen to make the sum (or 1-norm) of equal to 1. The distribution of initial values for uses that is random, initially uniformly distributed between 1 and 2, and then normalized to unit 1-norm. Also shown is the fit at different instants of time. These results demonstrate that provided the initial infected population distribution is small, the subsequent evolution obeys a simple one-dimensional description with a scalar variable f(t)

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

Continuum model

We consider a continuum limit of Eqs. (15) and (16), as . We replace the summation in Eq. (9) with integrals and assume and to be functions of and t. In the continuum limit, Eqs. (15) and (16) becomewhereIn 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 beand the fraction of the presently infectious population is understood to beAt the start of the pandemic, we have initial conditions that satisfyandby which we mean a tiny initiation of infection.

Introduction of a tiny infected population can, with suitable initial distributions of S, make the pandemic grow. A society that has strong social distancing will have decaying to zero rapidly with increasing and may not see strong growth of the infection. Conversely, a society where the population has greater social mixing, i.e., 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 and are nonnegative. It is in principle possible for 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 and turns out to be satisfied.

Motivated by the foregoing observations from numerical simulations, we now look for solutions of the formwhere can be determined from initial conditions. Substituting into Eq. (20) yieldsThe integral on the right hand side is purely a function of time t: let us call it g(t). Then, we haveAlso, Eq. (21) becomesWe multiply both sides above by and integrate with respect to . Definingand assuming the integral and the derivative on the left hand side can be interchanged, Eq. (25) becomesAny 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 from Eq. (27), we obtain:Separately, multiplying both sides Eq. (26) by , we obtainwhere we notice that the right hand side is integrable with respect to time, givingwhere 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 writeif it suits our purpose. Equation (28) becomes, upon one integration,where is a constant of integration and can be evaluated by setting when (at the time of infinitesimal initiation of infection) and is given by

For many initial distribution functions , G(f(t)) and can be evaluated in closed form. For example, for the uniform distributionwe havewhereandIn 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 infinite-dimensional 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 obtained from Eq. (33), and the finite network solution [Eqs. (15) and (16)] fitted to for different parameter values (see figure caption for details). The results match closely. From Fig. 5, it is evident that, if the initial conditions for the infected population are all sufficiently small, the subsequent evolution does not depend very much on them. This is because the interaction term smoothes out all those initial variations, and the dimensional collapse occurs, which is the main point of the paper. In particular, see the simulation results in Fig. 5a and b. In one case, the infective population is taken to have a smooth variation with respect to , and in the other case, it is a random variation with respect to . The results are essentially the same, and the match is nearly perfect.

Fig. 5

Comparison between the continuum solution of Eq. (33) and the discrete network solution [Eqs. (15) and (16)] with fitted using . a and . b and randomized as in Fig. 4. The history function for Eq. (33) is taken as , with is adjusted to match the results from [Eqs. (15) and (16)]. Other numerical parameters are mentioned in text boxes within subplots

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 .

Equation (31), written as Eq. (32), becomesNow, if we interpret as the probability density function of the random variable in the population and think of as the probability density function of the transformed variable in the same population, then it must be true thatwhich giveswhere 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 [29] of the distributed quantity .

If we expand in a series for small f, we obtainwhere the coefficients are moments of u, given byassuming these moments exist. In other words, is the population mean of , is the population mean of , is the population mean of , 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 haveIn Eq. (42), if we are interested in an outbreak that starts from an infinitesimal infection, such that when , then the constantconsistent with Eq. (34), leavingLinearizing Eq. (43) for small f, we obtainSubstituting , we findThis characteristic equation matches one studied in [20, 26], and the stability condition is known to be [recall Eq. (22)]In fact, since 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 [20, 26]. 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 (39) yieldsEquation (33) becomes (inserting to match initial conditions of when )Setting , we obtainwhich exactly matches Eq. (15) of [26].

Further, even for variable , the correpondence between Eqs. (33) and (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 , and (ii) is the same, up to a linear scaling, as the dynamics with a single fixed as studied in [26]. For larger outbreaks, higher moments of the distribution of begin to play a role, and the match with Eq. (49) deteriorates.

Finally, we close with an example where the required moments in Eq. (43) do not all exist. Let the probability density function of befor which, in Eq. (43), and . This distribution has a long tail. The integral G(f) for [recall Eq. (40)] can be found (using Maple) in terms of special functions, and the first few terms of an asymptotic series for small is as shown below:where is Euler’s constant, also called the Euler–Mascheroni constant. The series has an term, reflecting the unbounded third moment of u. Equation (51) can be used in a numerical solution of the DDE in Eq. (33) with Eq. (34), provided the outbreak is not too big. For comparison, straightforward solution of the governing equations can be tried with direct numerical integration of Eq. (33) combined with Eq. (39) truncated at some suitably large value of u (we used ). Some such solutions are shown in Fig. 6. The match is essentially perfect, which indicates that finiteness of moments of u, and hence the moments of are not needed for the underlying DDE to have useful solutions.

Fig. 6

An example with unbounded moment (see Eq. (50)). a on the horizontal axis governs stability. Each circle, labeled “numerical,” represents the saturation value observed in a numerical solution based on integration of Eq. (33) combined with Eq. (39). At every time step, the integral in Eq. (39) is evaluated numerically for the truncated domain . The continuous line labeled “analytical” is from Eq. (33) with Eq. (51) used directly. b Comparison of time responses , with “numerical” and “analytical” implying the same as in (a). The initial function used was

Discussion and conclusions

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 and 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 . 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 [30] 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 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 finite-moments 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 , or of and , 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.