A unifying nonlinear probabilistic epidemic model in space and time

Covid-19 epidemic dramatically relaunched the importance of mathematical modelling in supporting governments decisions to slow down the disease propagation. On the other hand, it remains a challenging task for mathematical modelling. The interplay between different models could be a key element in the modelling strategies. Here we propose a continuous space-time non-linear probabilistic model from which we can derive many of the existing models both deterministic and stochastic as for example SI, SIR, SIR stochastic, continuous-time stochastic models, discrete stochastic models, Fisher–Kolmogorov model. A partial analogy with the statistical interpretation of quantum mechanics provides an interpretation of the model. Epidemic forecasting is one of its possible applications; in principle, the model can be used in order to locate those regions of space where the infection probability is going to increase. The connection between non-linear probabilistic and non-linear deterministic models is analyzed. In particular, it is shown that the Fisher–Kolmogorov equation is connected to linear probabilistic models. On the other hand, a generalized version of the Fisher–Kolmogorov equation is derived from the non-linear probabilistic model and is shown to be characterized by a non-homogeneous time-dependent diffusion coefficient (anomalous diffusion) which encodes information about the non-linearity of the probabilistic model.


Results
The probabilistic model. In order to illustrate the rationale of the model we can look at Fig. 1 which shows the spread of the COVID-19 on a global scale from the 22th of February to the 20th of May 2020. A logarithm density plot of the number of infected people using a geo-referenced map is shown in Fig. 1 which highlights the east-west direction of the virus propagation from China to Europe and then to America, with an increasing propagation of the outbreak. Let be the surface of the earth we focus on. A point in the surface is denoted by x = (x, y) ∈ � . We assume there are positive definite functions ψ I t (x), ψ R t (x), ψ S t (x), ∈ L 1 (�, d 2 x) representing sub-probability densities and describing the state of the epidemic at time t. In particular, ψ I t (x) is the sub-probability density for the infected, ψ R t (x) the sub-probability density for the recovered (infected which turn into immune) and ψ S t (x) the sub-probability density for the susceptible. The probability density functions are assumed to be derivable with respect to the time variable. Moreover, we assume ψ S 0 (x) + ψ I 0 (x) = ψ S t (x) + ψ I t (x) + ψ R t (x) for every t and � (ψ S 0 (x) + ψ I 0 (x)) dx = 1 . As a consequence we are not considering births and deaths. It is worth noticing that the existence of ψ I t , ψ R t and ψ S t is a theoretical assumption; the actual infected distribution is just an instantiation of the random process associated to the infected distribution. Anyway, an approximation of ψ I t , ψ R t and ψ S t can be obtained through density estimation methods as for example the kernel density estimation 38 (KDE) which provides a differentiable probability density function.
We suggest a partial analogy with the interpretation of quantum mechanics 39 where the state of a particle is described by a probability cloud which is given by the square of the modulo of the wave function . In particular, |�(t, x)| 2 ∈ L 2 (R) is a probability density function and � |�(t, x)| 2 dx is interpreted as the probability that the particle is found to be in at time t. No further specification is required. Analogously, the state of the epidemic is defined by three probability clouds (ψ I , ψ S , ψ R ) and � ψ I t (x) d 2 x is interpreted as the probability that an infected is in the region ⊂ at time t. No further specification will be required in the present approach. The statistical interpretation of quantum mechanics 39 introduces the concept of ensemble of identical copies of the particle in order to explain the meaning of the function . To be more precise, describes an infinite ensemble of copies of the particle; each copy being prepared by the same experimental arrangement. The function | | 2 then gives the statistical distribution of the position of the particle among the copies of the conceptual ensemble. We can adapt such an interpretation in the present context by interpreting ψ I as follows: suppose there is a mechanism which encodes the nature of the epidemic as well as the contingent circumstances it is determined from; suppose such a mechanism can be used to produce an individual; he/she will be an infected with probability www.nature.com/scientificreports/ and � ψ I t (x d 2 x is interpreted as the probability that a copy of the ensemble is made of an infected located in the set . That provides a possible interpretation of the state of the epidemic as we have defined it and is helpful, we think, in order to explains the role that ψ I plays in the present paper. In quantum mechanics, the dynamics of the system is determined by a deterministic evolution equations for (the Schroedinger equation). Pushing further the analogy with quantum mechanics, the dynamics of the epidemic is given by the evolution equations for ψ I and ψ S [see Eq. (6) and (8)]. Note that a deterministic evolution for the probability density does not mean that the epidemic evolves deterministically. It is a probabilistic process.
One could ask what is the connection between the state of the epidemic as we have defined it and the actual observation of the epidemic. Suppose we are able to generate an individual as many times as we want as we discussed previously. Suppose we use such a generation procedure to provide a family of ensembles E 1 , . . . E k each one made of N copies of the individual according to the state ( ψ I t , ψ S t , ψ R t ). The ensembles E i can be considered as instantiations of the epidemic with N individuals. The distribution of the infected, σ I j (number of infected over surface), corresponding to the instantiation E j is in general different from the distribution of the infected σ I i corresponding to E i . Despite the number of infected (susceptible) in a given region will be different in any instantiation, the expected number of infected (susceptible) in a given region is meaningful from the statistical viewpoint and is given by . By the law of large numbers, the probability that the number of infected in E i differs from N � ψ I t (x) d 2 x can be made arbitrarily small by choosing N sufficiently large.
We remark that in the present paper the focus is on the probability density functions ψ I t , ψ S t , ψ R t which, by analogy with quantum mechanics, define the state of the epidemic. The approach is then phenomenological as opposed to mechanistic; we assume the existence of probability density functions describing the probabilistic distribution of infected, susceptible, and recovered (the state of the epidemic) and limit ourselves to describe their time evolution for which we postulate a deterministic law (see below). The main aim of the present paper is to show that this approach defines a non-linear probabilistic model from which many of the existing models can be derived providing at the same time a unifying framework. Moreover, the model could be used to describe the evolution of the infected spatial probability density in order to locate those regions of space where the probability of infection is going to increase or to decrease. That could be very helpful to governments as we remarked in the introduction and will be the topic of a future paper. Now, we proceed to define the evolution equation for the state of the epidemic. We provide a non-linear evolution for the sub-probability distributions.
In the following, α(t) denotes the probability rate that infected turn to immune; it is assumed to be continuous and such that, t 0 α(τ ) dτ = p(t) ≤ 1 , for all t ∈ R + , where p(t) is the probability that an infected turns into immune during the time interval [0, t]. The probability rate α can be used to connect ψ R and ψ I : The space dependence of α can also be included if necessary. Moreover we assume there are no births and deaths and no displacements of the individuals (see the comments after Eq. (3) and before 3) has been used to estimate ψ I from the data about the space distribution of the infected. The software uses the kernel density estimation method 38 to estimate ψ I . Then, the same software has been used in order to generate a sampling, i.e., a set of points on the surface distributed according to ψ I . That has been done at several instant of time. The image illustrates qualitatively the expansion of the infected population. The infected cloud on the ocean is an artifact. www.nature.com/scientificreports/ Eq. (5) for further details). The epidemic propagation is encoded into a kernel W(t, x 2 , x 1 ) which describes the transition probability rate of the infection from x 1 to x 2 at time t. In particular, is interpreted as the conditional probability that a susceptible in the neighbourhood centered in x 2 with radius is interpreted as the probability that a susceptible in the neighbourhood centered in x 2 with radius d 2 x 2 gets infected by the infected contained in the neighbourhood centered in gives the probability of new infected in the neighbourhood centered in x 2 with radius d 2 x 2 during the time interval [0, t]. Since in the same time interval there is a probability that the infected turn into recovered, we have The previous probabilities can be interpreted by resorting to the idea of conceptual ensemble we introduced above. In particular, one can think of the ensemble of infinite copies of a single individual with each copy generated according to the state of the epidemic by the procedure described above; each individual being fixed in his/her position with the susceptible that can become infected and the infected that can become recovered according to the probabilities introduced above. In the case of an instantiation of the epidemic with N individuals fixed in their positions and supposing N sufficiently large, the evolution of the population densities ( σ I , σ S , σ R ) can be approximated by using the previous probabilities (see also Fig. 3 on this last point).
Note that, at variance with the second member of Eq. (3), the first member does not depend on α . The role of α is to give the probability that infected turn into recovered while the first member gives the probability of new infected which depends on α only indirectly: a higher α means a lower ψ I τ . The flux diagram in Fig. 2 illustrates the variation of the sub-probability distributions during an infinitesimal time interval dτ.
By (3) we obtain, with, One can think of as the transition probability that a susceptible in 2 gets infected during the time interval (τ , τ + dτ ) if in x 1 there is an infected. Note that µ (·) (τ , x 1 ) is absolutely continuous with respect to (·) ψ S τ (x 2 ) dx 2 and W(τ , x 2 , x 1 ) is the the Radon-Nikodym derivative of µ (·) (τ , x 1 ) with respect to (·) ψ S τ (x 2 ) dx 2 . That seems quite natural since Evolution of the sub-probability densities during an infinitesimal time interval dτ . Here, www.nature.com/scientificreports/ in terms of transition probability and explains the role that the susceptible probability density ψ S plays in the definition of W.
It is worth remarking that the kernel is intended to provide a phenomenological description of the evolution. We are assuming there are no displacement of the individuals. We just assume the existence of the transition probability (1).
Supposing ǫ sufficiently small, we can use (4) to derive the evolution from the state at time t to the state at time t + ǫ, Then, we can derive an integro-differential equation for the infected sub-probability density ψ I t . Indeed In order to derive the differential equation for the susceptible sub-probability density we write the susceptible sub-probability density at time t + ǫ as a function of the susceptible sub-probability density at time t, where the first term in the second member is the sub-probability density at time t and the second term is the probability density of the new infected produced in x 2 in the time interval [t, t + ǫ] . Then, by using the Taylor expansion for ψ S t+ǫ (x 2 ) , we obtain It is interesting to note that an equation similar to Eq. (6) and describing the evolution of the total size of an infective population I(t, ξ) in a deterministic framework can be found in Diekmann 18 (see exercise 8.21, page 219). It is derived in a different framework where x) dτ and i is the expected number of new infected per unit of time. The kernel A(t, x, y) defines the evolution of i. In particular, where S denotes the density of susceptible and a particular form of the kernel is assumed, A(t, x, y) = β(x, y)e −ατ . Note that in the present paper i assumes the form Figure 3. The upper part of the figure illustrates a hypothetical evolution of the normalized density of infected ρ I = σ I /N (density of infected over the total number of individuals) from time t to time t + t and shows that during the time interval t it can evolve in different possible ways because of the probabilistic nature of the system. That is not the case of the sub-probability density ψ I whose evolution is unique (lower part of the figure). Note also that N � ψ I t (x) dx can be very different from the number of infected contained in at time t (compare, for example, the yellow curve in the upper part of the figure with the red curve in the lower part). Only in the case N is very large (deterministic limit of the model) the evolution of the normalized density ρ I follows, with a high probability, the evolution of the probability density ψ I and the probability that N � ψ I t (x) dx coincides with the number of infected in the region at time t is high. Going back to the upper part of the figure, it is worth remarking that the deterministic models choose one of the possible evolutions for the normalized density of infected, ρ I . For example, in the Kendall model 40 a system of differential equations for ρ I t is postulated; it is equivalent (up to some restrictions on the structure of the kernel: the kernel is not time dependent and depends on the difference x − y ) to the system of differential equations for the probability density ψ I t we derived in the present work; such a postulate is sound only if N is sufficiently large. www.nature.com/scientificreports/ � ψ S τ (x)W(τ , x, y)ψ I τ (y) d 2 y . An equivalent, less general, system of differential equations for the evolution of the normalized densities ρ I t = σ I /N and ρ S t = σ S /N of infected and the susceptible (density over the total number of individuals) in a deterministic model can be found in Kendall 40 . For N very large we recover a generalized version of the Kendall model since ρ I , ρ I , ρ R provide an approximation of the sub-probability densities ψ I , ψ S and ψ R . Therefore, Kendall's model can be considered as a particular deterministic limit of the probabilistic model we suggest (see Fig. 3 for more details). An equation similar to Eq. (8) can be found in the deterministic model proposed by Thieme 41 which is based on a different approach and where the kernel k differs from W and assumes a different meaning. In particular, i(t, x) = � t 0 ∞ 0 I(t − τ , c, y)k(x; dy dc dτ ) where i is defined as the infective influence, I is the density of the infected and c denotes the age class. Then, the differential equation for the susceptible takes the form x) which compared with Eq. (8) shows the differences between k and W. We remark that, at variance with the space models just cited, in the present paper, Eqs. (6) and (8) refer to a probabilistic model and take on a particular meaning which is connected to the evolution of the sub-probability density functions describing the epidemic. They are assumed to grasp the law governing the evolution of the epidemic and then to provide the exact description of its evolution. We are assuming the existence of such a law which is encoded into the kernel W. Moreover, the integro-differential Eqs. (6) and (8) are derived in a probabilistic framework based on conditional probabilities. The probabilistic nature of the model and the differences with deterministic models is depicted in Fig. 3 below.
Extending further the analogy with quantum mechanics we can think of W as the analogous of the Hamiltonian operator in quantum mechanics which describes the interaction between the different particles of the system. Here W takes into account the interaction between infected and susceptible.
Note that Eqs. (6) and (8) define a non linear system of differential equations. The probabilistic framework we introduced is very general and is shown to contain several well known models that can be derived from it. We start with a generalized version of the Fisher-Kolmogorov model (Eq. (14)). Next we discretize the space coordinate and obtain a model which is continuous in time and discrete in space (Eqs. (16) and (18)). By neglecting the space effects, we derive the SIR and SI models (Eqs. (29) and (31)). By discretization of the time coordinate we obtain a general time-discrete stochastic model (Eq. (20)) from which the SIR stochastic model can be derived (Eq. 22).

Derivation of a generalized version of the Fisher-Kolmogorov model. The Fisher-Kolmogorov equation
where ρ t (x) denotes the density (or frequency) of a population is a particular kind of drift-diffusion equation that has been used in order to describe the frequency of a mutant gene in a population 42 (Fisher 1937). Later, a more general version has been studied by Kolmogorov 43 et al. (1937). We derive a generalization of the Fisher-Kolmogorov equation for the infected probability density ψ I that, assuming the stochastic fluctuations can be neglected, can be used to derive the evolution of the infected density of a population of N individuals. Then, under some stronger hypothesis, we derive the 2-dimensional version of Eq. (9).
Let us suppose that W(t, x 2 , x 1 ) depends only on the distance between x 2 and x 1 and = R 2 . Then Eq. (5) becomes where z = x 2 − x 1 . Furthermore suppose that in the time interval ǫ , ǫW is different from zero only for very small distances, |z| , and consider the second order Taylor expansion of where H denotes the Hessian matrix. By replacing ψ I t (x 2 − z) in Eq. (10) by its Taylor expansion and considering that W is an even function of the components of the vector z = (z 1 , z 2 ) , we obtain where and (9)

once replaced in the previous equation, gives
By renaming the terms we obtain where η := (aψ S 0 (x 2 ) − α) and µ = a . Equation (14) is a two dimensional drift-diffusion equation which generalizes the Fisher-Kolmogorov equation and it is worth remarking that it has been derived from a non-linear probabilistic model.
A comment is in order since Eq. (14) describes the evolution of a probability density while the Fisher-Kolmogorov equation refers to the deterministic evolution of the population density. Nevertheless, assuming that the infected population is sufficiently large that the stochastic effects can be neglected, an analogous equation can be derived for the infected population density which, in this case, can be approximated by the expected number of infected ρ I t (x) = Nψ I t (x) . One obtains from (14), where, γ = µ Nη . We obtain the Fisher-Kolmogorov equation by assuming α = 0 (which is equivalent to neglecting the removed) and by assuming ψ S t (x 2 ) ( ρ S t (x) ) constant and uniform. Those are reasonable assumptions if one is interested to study a single patch for sufficiently small time-intervals during which the susceptible are assumed to be uniformly distributed and the infected can be neglected with respect to the susceptible which can therefore be considered constant in time (the probabilistic model becomes linear under these hypothesis).
Thus, the Fisher-Kolmogorov equation is derived from a linear probabilistic model (the probability density of the susceptible is constant) and this extends previous results obtained by Mollison 34 who showed that although the Fisher-Kolmogorov model is nonlinear, it is insufficient to grasp the full nonlinearity of the phenomena. In our opinion, it is relevant that a generalized version of the Fisher-Kolmogorov equation which is characterized by the time dependence of the diffusion coefficient (see Eqs. (14) and 15) is instead derived from a non-linear probabilistic model. That indeed suggests that the information about the non-linearity of the epidemic evolution could be fully encrypted in the time dependence of the diffusion coefficient D t (x) (anomalous diffusion); the latter being necessary for the probabilistic model to be non-linear (see Eq. (13)). It is worth remarking that anomalous diffusion has been recognized to be very relevant in many biological and physical systems. For example, it has been shown to be characteristic of the motion of single messenger RNA molecule in a living Escherichia coli bacterium 48 . A possible explanation for anomalous diffusion of the kind D ∝ t α in biological systems can be given by continuous-time random-walk models. The derivation of Eq. (14) from the non-linear probabilistic model seems to suggest another possible reason for anomalous diffusion, i.e., non-linearity of the probabilistic model. Moreover, Eqs. (13) seems to suggest that more general kinds of anomalous diffusion are possible. Those are problems deserving further investigation.  (17) is obtained by multiplying the discrete version of Eq. (7) by σ j . By using again the Taylor expansion of p I j (t + ǫ) and p S j (t + ǫ) and assuming ǫ sufficiently small, we obtain Space-time discretization and derivation of the SIR stochastic model. It is worth remarking that the discrete time version of Eqs. (17) provides a discrete time stochastic process. By choosing ǫ sufficiently small, is interpreted as the probability that a new infected is produced during the time interval ǫ . The process in now discrete in time with the time intervals that are multiples of ǫ.
Once we have discretized the time, the SIR stochastic model [44][45][46] (also known as the general stochastic epidemic) can be derived. Indeed, if for every region, we neglect the contribution from the other regions, W t (j, i) = δ ij β i (t) , we obtain Now, let us denote by n j and s j the number of infected and susceptible in the region j respectively and assume ǫ sufficiently small to have n j (t + ǫ) − n j (t) ∈ {−1, 0, 1} . Then, we can give the following interpretation of the probabilistic model: for every fixed j, n j and s j are random variables with expected values Np I j and Np S j respectively, the transition probability from the state (n j (t), s j (t)) to the state (n j (t) + 1, s j (t) − 1) is given by p S j (t) ǫβ j (t)p I j (t) while the transition probability from the state (n j (t), s j (t)) to the state (n j (t) − 1, s j (t)) is given by ǫαp I j (t) and the transition probability from the state (n j (t), s j (t)) to the state (n j (t), s j (t)) is given by 1 − p S j (t)ǫβ j p I j (t) − ǫαp I j (t) . Note that the transition probabilities are derived from the general model; they are not assumed by definition.
Derivation of SIR and SI models. Finally, we note that deterministic SIR and SI models can be obtained by neglecting the space effects. In order to show that, we assume N sufficiently large for the stochastic effects to be neglected and re-write Eqs. (18) in terms of the number of infected n j (t) = Np I j (t) and the number of susceptible s j (t) = Np S j (t) in the j-th region. We obtain, www.nature.com/scientificreports/ where, � t (j, i) = 1 N W t (j, i). If we neglect the contribution of the other regions, i.e., � t (j, i) ≈ δ ji c j (t) , (which means that we are neglecting the space effects on the evolution of the epidemic), the system of differential equations becomes which coincides with the equations for infected and susceptible of an SIR model with time-dependent parameters c j (t) , α(t) and with birth rate and death rate equal to zero (we assume that deaths and births compensate each other). SIR models are particular cases of the stochastic model we have introduced. If moreover, we neglect the removed, we obtain an SI model In particular, the relation between infected and susceptible is s j (t) = s j (0) − n j (t) where s j (0) is the number of susceptible at time t = 0 and the equation for the infected becomes where k = c j (t)s j (0) and a = s j (0) , which is a logistic equation with time-dependent parameters. Assuming c j (t) time-independent we obtain the Verhulst logistic Eq. 47 (Verhulst 1838) whose solution is the sigmoid function where q = a−n j (0) n j (0) and n j (0) is the number of infected at time t = 0.

Discussion
We have proposed a general space-time continuous probabilistic model where the state of the epidemic is given by three sub-probability density functions. The latter are interpreted by analogy with the statistical interpretation of quantum mechanics. Then, we showed that many important stochastic and deterministic models can be derived as particular cases. That should be helpful from the theoretical viewpoint since it can be used to analyze further the relationships between different modelling approaches as well as their limits. For example, it has been shown that the non-linearity of the Fisher-Kolmogorov equation is not sufficient to characterize the non-linearity of the stochastic phenomena it aims to approximate 33,34 . We strengthen and generalize such a result but also show that a general version of the equation characterized by anomalous diffusion (the diffusion coefficient is spacetime dependent) is more deeply connected to the non-linearity of the probabilistic model it is derived from. That is of general interest from the physical viewpoint as well as for biological applications since it concerns the relationships between non-linear stochastic and non-linear deterministic models. Moreover, anomalous diffusion plays a relevant role in the diffusion of single messenger RNA molecules in living cells and in many other physical and biological systems [48][49][50][51] . A possible explanation for anomalous diffusion of the kind D ∝ t α in biological systems can be given by continuous-time random-walk models. The derivation of Eq. (14) from the non-linear probabilistic model seems to suggest another possible reason for anomalous diffusion, i.e., non-linearity of the probabilistic model; it furthermore suggests that more general kinds of anomalous diffusion could be possible. A future work will be devoted to push forward such investigations.
Possible applications to Covid-19 epidemic concern the time evolution of the probability density ψ I t . The density ψ I t can be approximated by the kernel density estimation method (KDE) and its evolution is determined by Eqs. (6) and (8). The knowledge of such evolution, which is encoded into the kernel W, would allow probabilistic epidemic forecasting; it would allow to locate those regions of space where the probability of the infection is (25) = −Nα(t)p I t (j) + Np S j (t) n i=1 1 N W t (j, i)Np I t (i) (26) = −αn j + s j n i=1 � t (j, i)n i (27) s j (t) = ds j (t) dt = N dp S j dt (28) = −s j n i=1 � t (j, i)n i (29) n j (t) = −αn j (t) + s j (t)c j (t)n j (t) (30) s j (t) = −s j (t)c j (t)n j (t) (31) n j (t) = s j (t)c j (t)n j (t) (32) s j (t) = −s j (t)c j (t)n j (t) n j (t) = s j (0)c j (t)n j (t) − c j (t)n 2 j (t) = kn j (t) − k a n 2 j (t) n j (t) = a 1 + qe −kt www.nature.com/scientificreports/ going to increase. All of that depends on our capacity to estimate the transition kernel W. It is a kind of inverse problem that requires insight into the real scenarios. Concerning this last point, some preliminary work 52 that could be helpful and some applications 53 of the kernel density estimation method based on ideas that could be relevant, already exist but much more effort is necessary in order to estimate W starting from data. More mathematical work is required as well. That will be the topic of a future work.