Inference of a Susceptible–Infectious stochastic model

1.   Introduction

The onset of large epidemics in recent decades, such as the Severe acute respiratory syndrome (SARS) epidemic, the avian influenza (^[1], ^[2]), Ebola (^[3]), and the Covid-19 pandemic (^[4], ^[5]), have meant that a lot of mathematical models to study various infectious diseases have been developed (see, for example, ^[6], ^[7], ^[8], ^[9], ^[10]). The main aim of such studies is to forecast the dynamics of the disease, by using also suitable inference techniques in order to adopt appropriate containment policies (see ^[11], ^[12], ^[13]).

Generally in epidemic models, the population is divided into compartments that constitute a partition of it, that is, each compartment constitutes a subset of the population disjoint from the others and the union of all the compartments returns the whole population. Among the models that contemplate this compartmentalized philosophy, the susceptible-infectious-recovered (SIR) type and its derivatives stand out.

The study of disease propagation through these models, both in their deterministic and stochastic versions, has experienced a great boom in recent decades, giving rise to an extensive literature. The range of models considered is very wide, taking into account different points of view. For example, in ^[14] and ^[15], Kalman filtering techniques are applied to estimate the states of a discrete nonlinear compartmental model. In the stochastic environment, Markovian models occupy a prominent place. Within them, if we consider models with discrete states, continuous-time Markov chains have been widely used (see, for instance, ^[16], ^[17], ^[18], ^[19], ^[20], ^[21]). As for continuous-time Markovian models with continuous state space, diffusion processes arise naturally by introducing random environments (through a multidimensional Wiener process) into the systems of ordinary differential equations governing deterministic models.

From the classical SIR model, a number of increasingly complex variants have emerged. Among the lines on which the evolution of the models has been based, we can highlight the following:

● The partitioning of the total population into an increasing number of compartments, each of which obeys a particular situation of the individuals. Thus, in addition to the usual susceptible (S), infected (I) and recovered (R) individuals, there are others such as those with passive immunity (M), the exposed and uninfected (E) and deceased (D), among others. This has led to an increased complexity and difficulty in treating compartment models (see, for instance, ^[22,23]).

● The inclusion of the mechanisms of contagion of the disease, the so-called incidence function. The most usual function of this type is the one that relates susceptible individuals to infected individuals according to the law of mass action, but considering a constant contagion rate and/or rational functions of both groups of individuals to regulate the contagion scheme (^[24,25]).

● The inclusion of terms describing the effect of vaccination of individuals (including the possibility of restrictions in the vaccination process), as well as the possible effects of cross-infection (see, for instance, ^[26], ^[27], ^[28], ^[29])

● The emergence of new diseases. The effect of COVID19 on the increase in the literature on this type of model cannot be denied. Indeed, the analysis of the effects of the pandemic has given rise to numerous publications focusing on two lines of action: on the one hand, the study of the evolution of the disease in specific locations using existing models (see ^[30]) and, on the other hand, the appearance of new models (see, for example, ^[31]).

However, it should be noted that a very high number of these publications focus their interest on the description of the model and its analysis from a theoretical point of view (existence, uniqueness, non-explosion of the solution, stability, equilibrium problems and extinction of the disease). However, especially in stochastic models, the aspects derived from the estimation of the models are not yet well developed, resorting to simulations and/or Monte Carlo type methods to illustrate the validity of the proposed models. The main reason for this derives from the difficulty to know the probability distribution governing the dynamic evolution of the model (the more complex the more complicated the model). Some approximation has been made on simpler models (such as the susceptible-infectious-susceptible (SIS)) by means of approximations derived from Euler-Maruyama type discretizations (see, for instance, ^[32]), although this requires certain conditions on the model components (e.g. the incidence function). The situation can become even more complex if the rates involved in the model are not considered constant but time-dependent.

This paper therefore falls along the line of inference in epidemic models. The aim is the inference of a SI type model, which is essentially a "simplified" version of an SIR type model. Precisely, in the SI model it is assumed that an individual can be in one of only two states, either susceptible ($S$) or infectious ($I$). Although this model is quite simple, it is adept at capturing several types of diseases in which individuals remain infected for life (e.g., brucellosis in domestic and wild populations, fox rabies). Even the most famous AIDS has been modeled in the literature using an SI type model (see, for example, ^[33]).

For $t \geq t_0$, we denote by $S(t)$ the number of susceptible individuals, by $I(t)$ the number of individuals infected and by $K$ the total population size, where $K = S(t) + I(t)$ is constant. In similar models, infected individuals are lifetime infectious. We point out that here the size of the population is assumed to be constant, i.e., the birth and death rates of both populations S and I are assumed to be negligible. This is a very strong assumption, but reasonable if one observes the phenomenon for a limited time.

A model with two states is described via ordinary differential equations in the deterministic dynamics. More general models can be built making use of stochastic approaches based on the birth and death processes or diffusion processes (see, for instance, ^[6], ^[34,35], and references therein); they are more realistic but more complicated to analyze. Such models try to forecast the spread of the disease in terms of the total number of infected people and the duration of an epidemic. Moreover, they allow us to estimate suitable epidemiological parameters, as the transmission rate of the disease measured via the basic reproduction number, i.e., the expected number of infected cases directly generated by one case in a population where all individuals are susceptible or infected.

Let $I(t)$ and $S(t)$ be the sizes of the infected and the susceptible populations, respectively. We consider the deterministic SI model described by the following

where the transmission intensity function $\lambda(t)$ is a positive, bounded and continuous function of $t$. It follows that the population dynamics of the infected $I(t)$ can be described by the Pearl-Verhust logistic growth differential equation:

where the transmission intensity function $\lambda(t)$ is a positive, bounded and continuous function of $t$. The solution of (1.1) is

with

We note that $\lim_{t\to\infty} I(t) = K$, so the whole population is destined to become infected, hence the parameter $K$ identifies the carrying capacity of the infected population. We point out that the susceptible-infectious epidemic model is an extreme case of the more general models including recovered population and can be obtained from them by assuming that the time required to reach an immunity situation is infinitely long (see, for instance, ^[36]). The state $K$ is not reachable in finite time, so it is interesting to consider the time $T^*_m$ required to reach a fixed threshold $m$. When the process is time-homogeneous, i.e. the function $\lambda(t)$ is constant, one can determine the time $T^*_m$ such that $I(T^*_m) = m$. In particular, from (1.2) one obtains

being $I(t_0)$ the initial size of the infected population and $m\in(0, K)$.

To generate a stochastic diffusion process for $I(t)$, several approaches can be used, generally, they introduce stochastic elements into the equation (1.1) or into solution (1.2). These two approaches lead to different processes described by stochastic differential equations of Itô or Stratonovich type, respectively (see ^[37]). In this work, the first of the two approaches has been chosen.

In ^[35], a time-inhomogeneous diffusion process to model the size of the infected population in a stochastic environment is provided by starting from (1.2). For such a process, the authors study the probability distribution and derive closed form results for the first-passage time problem through a constant boundary to obtain the stochastic counterpart of the parameter $T^*_m$. However, such important issues as that of model inference were not addressed in this study.

The aim of the present paper is to provide the inference on the stochastic diffusion process built from (1.2). Indeed, such issue becomes fundamental to study the dynamics of infectious diseases and to measure their aggressiveness so as to make predictions about future infections.

The paper is organized as follows. Section 2 provides the stochastic model and its probability distribution, i.e. the transition probability density function (pdf) and the related moments. Section 3 contains the inference procedure based on the Generalized Method of Moments (GMM) for fitting the parameters and the unknown functions of the model. In Section 4, some simulation experiments are performed to validate the proposed procedure. An application to real data is also considered in Section 5. Some conclusions close the paper.

2.   The model

Under the assumption of random environment, let $\{X(t), t\geq t_0\}$ be the stochastic process describing the size of the infected population at time $t$, and we interpret $\Lambda(t|t_0)$ as the mean of a time-inhomogeneous Wiener process $\{Z(t), t\geq t_0\}$, described by the stochastic equation

where $W(t)$ is the standard Wiener process and

with $\sigma(t)$ being a positive, bounded and continuous function of $t$ describing the breadth of the random oscillations. Hence, Eq. (1.2) is generalized by the following stochastic equation:

As proved in ^[35], $\{X(t), t\geq t_0\}$ is a time-inhomogeneous diffusion process defined in $(0, K)$ described by the following stochastic differential equation:

where

are the infinitesimal drift and the infinitesimal variance of $X(t)$, respectively, and $x_0\in (0, K)$ is the initial size of the infected population.

As shown in ^[35], the process $X(t)$ can be transformed into a time-inhomogeneous Wiener process $Y(t)$ with drift and infinitesimal variance

by using the transformation:

We point out that $Y(t)$ is a Gauss Markov process with mean

and covariance function

Further, since the transition pdf $f_Y(y, t|y_0, t_0)$ of $Y(t)$ is normal with mean $y_0+\Lambda(t|t_0)$ and variance $V(t|t_0)$, we can obtain the transition pdf of the process $X(t)$:

Moreover, the transition distribution function of $X(t)$ is

where ${\rm Erf}(z) = \frac{2}{\sqrt \pi}\int_0^z e^{-u^2}du$ is the error function.

The conditional median $\mu[X(t) |X(t_0) = x_0]$ of the process $X(t)$ can be obtained from (2.11) by imposing that $F_X(x, t|x_0, t_0) = 1/2$ for $t\geq t_0$, soone has:

Moreover, for $m = 1, 2, \ldots$, the $m$-th conditional moment of $X(t)$ is:

Such analytical results will be used in the following section to make inference for the process $X(t)$, to provide a fitting procedure based on GMM.

3.   Inference

The transformation (2.7) is the basis of our estimation procedure for the functions $\lambda(t)$ and $\sigma^2(t)$. The idea is to estimate such functions by making inference on the transformed Wiener process $Y(t)$ in order to fit $\lambda(t)$ and $\sigma^2(t)$. In particular, from (2.8) and (2.9) we obtain:

In the following, we assume that the carrying capacity $K$ is known since it represents the total size of the population, whereas the functions to be estimated are $\lambda(t)$ and $\sigma^2(t)$, for $t$ belongs to the observation interval $[t_0, T]$.

We consider a discrete sampling of the process (2.4) based on $d$ sample paths observed at the times $t_{j}$, with $j = 1, \ldots, n$. For $i = 1, \ldots, d$, let $x_{ij}$ be the observed values of the $i$-th path at times $t_{j}$, and the values $x_{i1}$ represent the initial point of the sample paths.

The inference procedure is illustrated in the following:

● From the observed data $x_{ij}$ for $i = 1, \ldots, d$ and $j = 1, \ldots, n$ of the process $X(t)$, obtain the points $y_{ij}$ as the following:

Such points can be considered as observations of the Wiener process $Y(t)$ given in (2.6).

● From the data $y_{ij}$ for $i = 1, \ldots, d$, obtain the sample mean $\mu_j$

and the sample covariance $\nu_j$ between two subsequent observations $Y(t_{j-1})$ and $Y(t_{j})$:

● Interpolate the values $\mu_j$ for $j = 1, \ldots, n$ and $\nu_j$ for $j = 2, \ldots, n$. Let $\widehat M(t)$ and $\widehat \Sigma(t)$ be the functions obtained.

● Evaluate the derivatives of $\widehat M(t)$ and $\widehat\Sigma(t)$.

● From Eq. (3.1) and (3.2), obtain the estimate of $\lambda(t)$ and $\sigma^2(t)$ as follows:

We point out that the procedure just illustrated combines a GMM estimation with the interpolation of the sample mean and covariance points, so the consistence of the estimators (3.6) of $\lambda(t)$ and $\sigma^2(t)$ derives from the consistence of the GMM estimator and from the uniform convergence of the interpolation method, for example, in our analysis we can consider cubic spline interpolation.

Finally, we note that in real applications, it could happen that the observed paths do not reach the carrying capacity since the phenomenon is observed before $K$ is achieved. In these cases, we argue that an a priori rough estimate of the parameter $K$ can be obtained by starting from the maximum observed point. Such an estimate can be used to obtain the points $y_{ij}$ for $i = 1, \ldots, d$ and $j = 1, \ldots, n$ from which, by using the previous procedure, we fit the functions $\lambda(t)$ and $\sigma^2(t)$. An improvement of the estimate of $K$ could be then obtained by looking at the conditional median function (2.12) and by implementing a step-by-step procedure similar to that proposed in ^[11]. This topic will be the subject of future investigations.

4.   Some simulation experiments

In order to evaluate the suggested procedure, in the following we consider several simulation experiments in which $d$ sample-paths of $X(t)$ are simulated, each including $n$ equally spaced observations in $[t_0, T]$ with $t_0 = 0$, $T = 50$, $t_i - t_{i-1} = \Delta = 0.01$ $(i = 1, \ldots, n)$ and $x_0 = 20$. The estimation of the unknown functions is replicated $N = 500$ times. In all the cases we choose the carrying capacity $K = 200$. In Section 4.1 we consider the case in which both the functions $\lambda(t)$ and $\sigma^2(t)$ are time independent, whereas in the Section 4.2 one or both of them are continuous functions of the time.

4.1.   The time homogenous case

Suppose that $\lambda(t)$ and $\sigma^2(t)$ in (2.5) are constant functions. In particular, in this simulation study, we consider

In this case, the process $X(t)$ in (2.5) is time homogeneous, so the inference for the parameters $\lambda$ and $\sigma^2$ can be made by means of the classical MLE. Indeed, it is easy to obtain the MLE for $\lambda$ and $\sigma^2$ by looking at the log-likelihood function of the process $Y(t)$ in (2.6):

In particular, we obtain the following estimator:

In Figure 1 we compare the results obtained via the MLE method with those obtained via our procedure in the following denoted by GMM. In particular, in Figure 1 on the top the box-plots of the 500 estimates of $\lambda$ (on the left) and of $\sigma^2$ (on the right) are shown. We can see that the estimates of $\lambda$ by means of MLE and GMM have quite the same distribution ranging from $0.46$ to $0.51$. Differently, the estimates of $\sigma^2$ via MLE present a very low variability compared to those of the GMM, resulting that MLE is to be preferred in this case. However, we point out the MLE assumes that the parameters $\lambda$ and $\sigma^2$ are constant and so it looks at the log-likelihood as a function of such parameters; instead, our procedure can be applied in the general case in which we do not have such information and the parameters generally depend on time. On the bottom of Figure 1 the kernel density of the MLE and GMM standardized estimates of $\lambda$ (on the left) and $\sigma^2$ (on the right) are shown. The red curves represent the density of MLE, whereas the black curves refer to GMM. We can observe that the sampling distributions of both the estimators and both the methods are quite symmetric and superimposable.

Finally, let

be the mean relative error (MRE) of the estimators $\widehat\lambda$ and $\widehat\sigma^2$. In Table 1 the MREs are shown for different choices of the parameters $\lambda$ and $\sigma^2$, i.e. $\lambda = 0.4, 0.6, 0.8, 1.2$ and $\sigma^2 = 0.05, 0.1$, and for the two considered methods. The MREs obtained via the MLE and via the GMM are generally comparable, especially for $\lambda$, even though the MLE method provides lower errors as expected by seeing the box plots in Figure 1. This one is due to the strong assumption of the constancy of the parameters. Moreover, the MREs increase as $\sigma^2$ increases for both the estimators and for both the methods.

4.2.   The time inhomogenous case

In this section we consider the general case in which one or both of the functions $\lambda(t)$ and $\sigma^2(t)$ are time dependent. In the following we analyze three different cases:

${\bf a.}$ $\lambda(t) = 0.4+\sin t, \qquad \sigma^2(t) = \sigma^2 = 0.1$;

${\bf b.}$ $\lambda(t) = 0.4+\sin t, \qquad \sigma^2(t) = 0.1+0.01(1-e^{-2 t})^2$;

${\bf c.}$ $\lambda(t) = \lambda = 0.4, \qquad\quad \sigma^2(t) = 0.01(1.2+\sin t).$

We note that in all the cases the chosen functions are continuous and bounded. In particular, we choose $\lambda(t)$ constant or periodic to consider a seasonality effect of the infectious disease. Further, in the Case a and Case b, $\lambda(t)$ periodically becomes negative, including situations in which a period of growth in the infection rate is followed by a period of regression of it. Three different choices are instead made for the function $\sigma^2(t)$: constant, asymptotically constant and periodic.

The results for the Case a are shown in Figure 2 for the functions $\lambda(t)$ (on the left) and for $\sigma^2(t)$ (on the right). The red curve is the true function and the blue curve is the mean of the 500 obtained estimates $\widehat\lambda_i(t)$ and $\widehat\sigma^2_i(t)$, i.e., $\widehat\lambda(t)$ and $\widehat\sigma^2(t)$. The black lines represent the observed confidence interval for the functions $\lambda(t)$ and $\sigma^2(t)$, respectively, obtained as

where

are the standard deviations of the estimators. We can observe on the left of Figure 2 that $\widehat \lambda(t)$ fits very well in the true function $\lambda(t)$, and also the amplitude of the confidence interval is small due to low values of the standard deviation. For the function $\sigma^2(t)$, we can see the the estimate $\widehat \sigma^2(t)$ decreases to the true value $0.1$ starting from the time $t = 10$ and the confidence interval always contains the true value.

In the Case b we consider the same periodic function of the Case a for $\lambda(t)$, whereas $\sigma^2(t)$ is a time dependent function that asymptotically tends to the constant value $0.11$ and this value is quickly reached. This is the reason for which the results for the Case b, illustrated in Figure 3 for the function $\lambda(t)$ (on the left) and for $\sigma^2(t)$ (on the right), seem similar to the results in the previous case (Figure 2), although both the functions depend on time.

Very interesting are the results for the Case c shown in Figure 4. Here, the constant function $\lambda(t) = 0.4$ is estimated, by using our procedure and by means of a function $\widehat\lambda(t)$ showing a periodicity. Clearly, such periodicity is due to the presence of $\sin t$ in the function $\sigma^2(t)$. In this direction, we have to point out that the proposed procedure for estimating $\lambda(t)$ and $\sigma^2(t)$ is nonparametric and the only assumption on the unknown functions is boundness of such functions. So, in the Case c, the procedure, looking at the moments, and in particular to the mean and the covariance functions, captures the periodicity in the model and associates it to both the functions $\lambda(t)$ and $\sigma^2(t)$. Anyway, the mean estimate $\widehat \lambda(t)$ is very close to the constant function $\lambda(t) = 0.4$. Also, the function $\widehat \sigma^2(t)$ fits very well in the true function $\sigma^2(t) = 0.01(1.2+\sin t)$ as shown in Figure 4 on the right, although the confidence bands are increasingly further apart as time increases. These observations open the way to the possibility of finding a tool to discriminate between different estimated models. From this perspective, informational divergence could be a criterion to be used in a future study.

5.   An application to real data

In this section, we apply our estimation procedure to the dataset twentymeas included in the R package tsiR (see ^[38]). It contains biweekly data (IP = 2) related to measles infection for twenty locations in England from 1944-1964 and was studied in ^[39]. We point out that the application presented here is primarily for illustrative purposes of how the proposed procedure can be used with real data. People infected by measles become immune so more complicated models, such as SIR and susceptible- infected- recovered-susceptible (SIRS), may be more adequate. However, even if simplified, we think that the SI model can also be used in this case. In fact, it is true that individuals who recover from measles become immune, but they remain infected, in a certain sense. The only difference is that they cannot infect other individuals. In our opinion, this aspect is taken into account with the fact that the transmission intensity function depends on time; in particular, we expect that $\lambda(t)$ becomes very small as time increases, showing that the disease is transmitted less frequently as recovered people have become immunized.

For our analysis we first consider the cumulative number of the infected for each location and we normalize this number by using the maximum number of population size, identified in the variable "pop" of the dataset. In Figure 5, the sample paths of the infected population and the normalized sample paths are shown. Here, we consider each normalized time series of the infected people as a sample path of a same diffusion process X(t) modeled via (2.3). From the normalized process data (on the right of Figure 5), we fix $K = 0.25$, i.e., the asymptotic infected population is $25\%$ of the total population. In Figure 6, the estimated function $\widehat \lambda(t)$ (on the top) and $\widehat\sigma^2(t)$ (on the bottom) are plotted for the whole period of observation on the left, and for the period 1946 to 1965 on the right. By looking at $\widehat \lambda(t)$, we can see that $\widehat \lambda(t)$ presents a sharp decrease in the first period of observation due to its high initial followed by a rapid tendency to assume constant behavior over time. Further, in this second period, the behavior of $\widehat \lambda(t)$ shows a certain seasonality, as expected in infectious diseases. Also on the bottom a higher value of the estimate $\widehat\sigma^2(t)$ is observed in the initial period, and after that it continues to decrease asymptotically to 0.

6.   Conclusion

We have considered a stochastic diffusion process for the time-inhomogeneous deterministic SI epidemic model and we have proposed an estimating procedure to inferring the considered process. A relevant issue concerning the model under consideration is the fact that it works quite well in situations where the total population size is large. In situations where the size is sufficiently small, several authors recommend the use of discrete state space processes. However, it should be noted that the proposed estimation procedure allows dealing when the population size is not too large. An example of this can be found in the example described at the beginning of Section 4, where the total size does not differ so much from the original (200 vs. 20).

Another interesting aspect is related to the assumption that the total population size is constant, which is usually a generalized assumption in this type of model. This is due to the fact that they reflect a study carried out over a period of time that is not too long to consider relevant variations in the size of the population. One line of future work is to consider a population growth pattern, either logistic, Gompertz, or other growth models.

To estimate the functions $\lambda(t)$ and $\sigma^2(t)$ that characterize the process, we have proposed a procedure based on the GMM method combined with the interpolation of the sample mean and covariance points, so the consistence of the estimators derives from the consistence of the GMM estimator and the uniform convergence of the used interpolation method. It should be noted that the proposed procedure does not make any assumptions about the functional form of the unknown functions. The only necessary assumption is that the functions involved are continuous and bounded in the observation interval. Several simulation studies have been carried out which demonstrate the validity of the proposed methodology.

The results for both the estimates obtained in the considered simulation studies seem very close to the "true" functions. In particular, for the time homogeneous case, we have compared the results obtained via the MLE procedure with those ones obtained with our procedure: concerning the parameter $\lambda$, representative of the infection rate, the estimates are very close with those obtained via our method. Some differences are instead found for the parameter $\sigma^2$, related to environmental variability, in which better performances are shown for the MLE due to the strong assumption of the constancy of the parameters. However, the MREs for the two methods and for the parameters are in all the cases comparable. Other simulation cases have been analyzed in which the functions $\lambda(t)$ and/or $\sigma^2(t)$ are time-dependent, with particular reference to situations in which seasonality effects of the dynamics of the infection are included. Finally, we have applied our estimation procedure to the dataset twentymeas included in the R package tsiR (see ^[38]). It contains biweekly data (IP = 2) related to measles infection for twenty locations in England from 1944-1964. Regarding the estimates obtained for $\lambda(t)$ and $\sigma^2(t)$ via our method, we observe a sharp decrease in the first period of observation due to their high initial values, followed by a rapid tendency to assume constant behavior over time. Further, in this second period, the behavior of the transmission intensity function shows a certain seasonality, as expected in infectious diseases. Instead, the estimate of $\widehat \sigma^2(t)$ rapidly decreases to a constant value without showing any periodicity.

Conflict of interest

The authors are special issue editors for Mathematical Biosciences and Engineering and were not involved in the editorial review or the decision to publish this article. All authors declare that there are no competing interests.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research is partially supported by MUR-PRIN 2022, project 2022XZSAFN "Anomalous Phenomena on Regular and Irregular Domains: Approximating Complexity for the Applied Sciences" (Italy), by MUR-PRIN 2022 PNRR, project P2022XSF5H "Stochastic Models in Biomathematics and Applications" (Italy), by PID2020-1187879GB-100 and CEX2020-001105-M grants, funded by MCIN/AEI/10.13039/501100011033 (Spain). G. Albano and V. Giorno are members of the GNCS-INdAM.