Positivity and Boundedness of Solutions for a Stochastic Seasonal Epidemiological Model for Respiratory Syncytial Virus (RSV)

In this paper we investigate the positivity and boundedness of the solution of a stochastic seasonal epidemic model for the respiratory syncytial virus (RSV ). The stochasticity in the model is due to fluctuating physical and social environments and is introduced by perturbing the transmission parameter of the seasonal disease. We show the existence and uniqueness of the positive solution of the stochastic seasonal epidemic model which is required in the modeling of populations since all populations must be positive from a biological point of view. In addition, the positivity and boundedness of solutions is important to other nonlinear models that arise in sciences and engineering. Numerical simulations of the stochastic model are performed using the Milstein numerical scheme and are included to support our analytic results. 1 Universidad Los Andes, gcarlos@ula.ve, Mérida, Venezuela. 2 Universidad de Córdoba, aarenas@correo.unicordoba.edu.co, Montería, Colombia. 3 Universidad EAFIT, mcogollo@eafit.edu.co, Medellín, Colombia. Universidad EAFIT 95| Positivity and boundedness of solutions for a stochastic seasonal epidemiological model for respiratory syncytial virus (RSV)


Introduction
Seasonal diseases in humans are inherent in the organic growth of man from infancy to old age. Some examples of such diseases are measles, diphtheria, chickenpox, cholera, rotavirus, respiratory syncytial virus (RSV ), vector-borne diseases including malaria and even sexually transmitted gonorrhoea [1]. In the modeling of the transmission of seasonal diseases, several nonlinear models of ordinary differential equations have been used [2], [3]. In these models, the variables commonly represent subpopulations of susceptibles (S), infected (I), recovered (R), latent (E), transmitted disease vectors, and so forth. Generally, the most important term is the nonlinear term βSI which comes from the law of mass action and where β is the transmission parameter of the disease [2]. In these models, periodic continuous |96 Ingeniería y Ciencia autores functions β(t) (called sometimes seasonally-forced functions) incorporate the seasonality of the spread of the disease in the environment [4], [5], [6], [7]. As an example of a seasonally-forced function, many authors use the expression β(t) = b 0 (1 + b 1 cos(2π(t + ϕ))), where b 0 > 0 is the baseline transmission parameter, 0 < b 1 ≤ 1 measures the amplitude of the seasonal variation in transmission, and 0 ≤ ϕ ≤ 1 is the normalized phase angle. Other examples of seasonal functions can be found in [5], [8], [6], where the models have a time dependent transmission rate. Recently, in [9], [10], [11], [12] deterministic epidemiological models have been studied, where the contact rate depends on more complex variables and flow parameters of one subpopulation to another depend on time. However, these models do not take into account the inherent randomness associated with the time variation of the data. Moreover, environmental fluctuations are some of the most important effects in real world systems. A large portion of natural phenomena does not follow a deterministic law exactly, but rather oscillates randomly around some average value [13].
In the seasonal epidemic models there are two environmental fluctuations that can be considered: the first is the natural and well recognized fluctuation of the seasonal transmission rate which is modeled by the seasonally-forced function β(t) and the second is the smaller random perturbations affecting the seasonal transmission rate which are modeled by Gaussian white noise. These last fluctuations arise from uncertainties and variations in environmental, demographic and all other parameters involved in the model system [14], [15], [16].
Thus, in this paper we propose to examine the effects of the stochastic process solution of the mathematical model presented in [17]. The stochasticity or noise is introduced on the baseline transmission parameter of the model given in [17], with the parameter perturbation technique, which is a standard technique in stochastic population modeling [14], [18], [19]. In addition it should be mentioned that several stochastic models have been applied to several issues including epidemics, for instance [20], [21], [13], [22], [23], [24]. Other types of stochastic models, such Markov's chain and graphs can be applied for several diseases with different characteristics [25]. Previous studies regarding positivity and boundedness of solutions of stochastic differential equation models have been developed without seasonal transmission parameters. In this paper we are interested in the introduction of stochasticity in the seasonal parameter in order to provide some additional degree of realism when compared with their deterministic counterparts [21], [26]. The main reason to introduce environmental stochasticity is that this can be a driving force that may change the deterministic dynamics of models [5], [27]. In addition, the positivity and boundedness of solutions are important to other nonlinear models that arise in sciences and engineering where uncertainty is included [28]. Thus a very important issue that arises when stochasticity is introduced in the seasonal epidemic model is positivity and boundedness of the solution.
For stochastic modeling different types of perturbations or randomness can be included in the model [25], [28]. However, in the art of modeling it is necessary to investigate which type of perturbation is more suitable for the particular disease. This issue is not a straightforward task and needs to be tackled with different statistical tools. However, each disease may have a different environmental stochasticity depending on the region. Historically, stochasticity has been introduced in many models of stock prices, economics, queueing theory, engineering, and biology. Nevertheless we do not rely on only this fact. For example, based on real hospitalization data of the respiratory syncytial virus (RSV ), we believe that the transmission of this type of disease is affected by several physical and social variables [29], [30]. In fact RSV is an illness for which the timing and severity of outbreaks in a community vary from year to year [6], [31], [29]. Here we assume that the transmission rate parameter of the RSV can be assumed equal to an average value plus a small time fluctuating term and this term follows a normal distribution with mean zero [32]. Many parameters have variations of this type. Possible explanations of this fact are a noisy environment due to several factors, such as temperature, humidity, pollution, transport, population and others. Our approach is to perturb the transmission rate parameter of the RSV by a Wiener process. Thus, from a theoretical point of view we |98 Ingeniería y Ciencia autores can rely on the existing framework for ordinary differential equations dealing with Gaussian white noise perturbations [33]. Additionally it is important to mention that other parameters of the model such as birth rate can be also perturbed stochastically but generally, in epidemic models, the most sensitive parameter is the transmission parameter of the disease [34], [32], [35]. However, future work may include the study of positivity and boundedness of solutions when all parameters are perturbed jointly.
This paper is organized as follows. Section 2 introduces a stochastic model for transmission of seasonal epidemic diseases in a population, where some stochasticity is introduced on the baseline transmission rate. In Section 3 we show the positivity of the solutions of the proposed stochastic model. Section 4 is devoted to the analysis of the existence and dynamical behavior of the solutions in the disease free state. The the property of stochastically ultimate boundedness is proved in Section 5. In Section 6, numerical simulations using some parameter values from clinical data of hospitalized individuals due to the respiratory syncytial virus RSV of the Spanish region of Valencia are presented in order to support our previous analytical results. Finally, in Section 7 discussion and conclusions are presented.

Derivation of the stochastic model
In [17] a generalized epidemic SIRS seasonal model dealing with population proportions of susceptibles S(t), infected I(t) and recovered R(t) has been presented and the authors showed the existence of periodic solutions, where the contact rate is a continuous periodic function β(t). The system is a first order ordinary differential equation of the formṠ where: 1. The population is divided into three classes: ing.cienc., vol. 13, no. 25, pp. 95-121, enero-junio. 2017.
• Proportion of susceptibles S(t), who are all individuals that have not the virus, • Proportion of infected I(t), being all the infected individuals having the virus and able to transmit the illness, • Proportion of recovered R(t), who are all the individuals not having the virus and with a temporary immunity.
2. The instantaneous birth rate µ(t) > 0 for all t ≥ 0, and equal to the instantaneous death rate, and we take it to be a continuous function.
3. For all the classes it is assumed the birth rate µ(t) is the same, because we assume that deaths associated with the disease are small.

The contact transmission rate (called seasonally-forced function)
is a function β(t) between classes S(t) and I(t), and is a continuous periodic function, which satisfies It is very common to approximate β(t) = b 0 (1+b 1 cos(2π(t+ϕ))), where b 0 > 0 is the baseline transmission parameter, 0 < b 1 ≤ 1 measures the amplitude of the seasonal variation in transmission and 0 ≤ ϕ ≤ 2π is the normalized phase angle [6].
5. The instantaneous per capita rate of leaving the infective class I(t) is called ν(t) and the instantaneous per capita rate of recovered class R(t) is called γ(t), which are non-negative, continuous and bounded functions in It is assumed that γ(t) < ν(t), for all t ≥ 0, i. e., the per capita rate of recovery is smaller than the per capita rate of leaving the infective and also that µ u , ν u , γ u , and µ l , ν l , γ l are positive real numbers defined by It should be mentioned that generally the recovery parameter ν(t) and the loss of immunity parameter γ(t) are assumed constant but here we include a more general case where these parameters can vary with time. On the other hand, birth and death rate µ(t) are assumed varying with time in other epidemic models where the size of the population is varying with time [2].
Fluctuations in the environmental conditions are some of the most important factors affecting real world systems. A large part of natural phenomena do not follow deterministic laws exactly, but rather oscillate randomly around some average value [13]. Therefore, in this paper we propose to examine the effects of the introduction of environmental noise on the positivity and boundedness of the stochastic process solution of the mathematical model proposed in [17]. The stochasticity or noise is introduced by perturbing the baseline transmission parameter b 0 in the model. This is a common technique used in stochastic population modeling [14], [18], [19], [36], [37] and it allows us to get analytical results [38], [18], [19], [39], [36], [37].
Let us now consider the model (1) with the perturbation on the baseline transmission parameter b 0 given by white noise. The use of Gaussian white noise is a good hypothesis in this model since it is assumed that the baseline transmission parameter oscillates randomly around some average value, due to some time varying disturbances. It is important to remark that this is the usual way to consider fluctuation in stochastic population dynamics [13]. In addition, the transmission parameter seems to be the parameter with faster variability since it depends on many variables such as temperature, humidity, pollution, transport and population that can vary every single day. On the other hand, the death and birth rates vary much more slowly and the recovery rate varies much less since is an intrinsic value of the disease.
In this way, if a stochastic perturbation is made on the baseline transmission parameter b 0 as in [20], [21], [23], an Itô type stochastic ing.cienc., vol. 13, no. 25, pp. 95-121, enero-junio. 2017. differential system of the following form is obtained: where is an Itô process, f is the continuous deterministic component or drift coefficient, g is the continuous random component or diffusion coefficient [33], defined by f : Thus, f is an d-vector valued function, g is an d×m matrix-valued function, and W (t) is a m-dimensional stochastic process having scalar Wiener process components. For our particular case, m = 1 and d = 3. Here, the perturbation to the baseline transmission parameter b 0 is introduced in the following form where α ∈ R is the intensity of the noise and η(t) is defined as the formal derivative of a standard Wiener process (a formal derivative because W (t) has probability one of being nondifferentiable). Since the Gaussian white noise has mean zero and the intensity of noise is small in comparison to the baseline transmission parameter b 0 , one gets positive values forb 0 . Since negative values make no sense, the very few negative values that may be produced during the simulation will be discarded. The new transmission rate is given by: Thus introducing the perturbation on equation (1) and multiplying by dt, one gets the following stochastic differential system given by with initial data S(0), I(0), R(0) T ∈ R 3 + .
Since (3) represents a population system, it is important that we do not obtain negative values. Thus, we must first prove the positivity of the solutions. Unless otherwise specified, let (Ω, F , P) be a complete probability space with a filtration {F t } t≥0 satisfying the usual conditions (i.e. it is increasing and right continuous while F 0 contains all P-null sets). Let W (t) be the one-dimensional Brownian motion defined on this probability space. Let X(t) = (x 1 (t), x 2 (t), x 3 (t)) T and The following theorem guarantees the existence and uniqueness of the solution of system (3). for t ≥ 0 and the solution will remain in R 3 ++ with probability 1.

Local behavior
In this section we analyze the existence and dynamic behavior of the solutions in the disease free state.

Existence of the equilibrium disease states
The following theorem shows the existence of a disease-free equilibrium point.
Theorem 4.1. If β u < ν l , then the system (3) is globally asymptotically stable, in the sense that for any initial value X 0 ∈ Ω, the solution X(t) will tend to the equilibrium point (1, 0, 0) asymptotically with probability 1.
Proof: We consider the second equation of system (3). Thus, for I(t), if I(0) > 0, we can use the function V 2 (t) = ln(I(t)), with I(t) ∈ (0, 1). Applying Itô's formula, one obtains that Then Simplifying and dividing the above inequality by t > 0, it follows that i.e., Using the fact that lim and from (10) we have that Therefore, if follows that lim

Behavior around of disease free
The following theorem, shows that the solution of system (3), is oscillatory around of point (1, 0, 0), if β u < ν l . Furthermore, the intensity of this oscillations depends of α.

Global behavior
Theorem 3.1 and Corollary 3.1, show that the solution of the system (3) with a positive initial value in Ω, will remain positive in Ω. The properties of positivity and nonexplosion are essential for a population system. Once they have been established, is very important to discuss some other properties of the solution of system (3). From a biological point of view, due to the limitation of resources, the property of stochastically ultimate boundedness is more desirable than the nonexplosion property. First, we show that the solution of system (3) is asymptotically bounded. Thus, we show first an important result, and then the stochastically ultimate boundedness follows directly.

Numerical simulations
In this section, numerical simulations using some parameters values from clinical data of hospitalized individuals due to the respiratory syncytial virus RSV of the Spanish region of Valencia are presented in order to support our previous theoretical results. In addition, we illustrate the effects of introducing stochasticity by means of numerical simulations. At first the deterministic model of RSV is simulated and later the stochastic model. The introduction of stochasticity into the deterministic model helps to understand the differences between the deterministic model and the real data of RSV . Thus, the stochastic model can be seen as a more sophisticated model. The virus RSV is the cause of acute respiratory infections in children younger than 2 years old, mainly bronchiolitis and pneumonia.
RSV is spread from respiratory secretions through close contact with infected persons or contact with contaminated surfaces or objects. This virus has been known since 1957, but only recently has the adult pathology been established. It is also the cause of 18% of the hospitalizations due to pneumonia in adults older than 65 [41]. RSV is a seasonal epidemic with minor variations each year and coincides in time with other infections as influenza or rotavirus, producing a high number of hospitalizations and, consequently, a saturation of Public Health Systems [42], [43], [44], [45]. RSV also causes repeated infections throughout life, usually associated with moderate-to-severe cold-like symptoms; however, severe lower respiratory tract disease can occur at any age, especially among the elderly or among those with compromised cardiac, pulmonary, or immune systems.
For the deterministic case model the parameter values have been taken from real clinical data of RSV hospitalizations of children younger than 4 years old during the period January 2001 to December 2004 from the CMBD (Basic Minimum Data Set) database of the Spanish region of Valencia [32], [44], [46]. This data was collected weekly taking into account other diseases such as bronchiolitis and pneumonia. In addition some unknown parameters were estimated using the least-squares fitting procedure and using the hospitalization data. The mean square error obtained is less than 0.002. The graphical representation of the model fitting can be seen in Figure 1.
As can be seen graphically, the deterministic epidemic model (1) agrees well with the observed epidemic data, but it is important to notice that the hospitalization data is not exactly periodic as the deterministic model predicts [17]. Thus, stochastic perturbations on parameters are a good explanation of these differences between the mathematical model and data, since without their use the simulation of the model will be exactly periodic. More details that include practical issues of the RSV mathematical modeling such as confidence intervals, parameter estimation, and numerical methods can be seen in [32], [42].
ing.cienc., vol. 13, no. 25, pp. 95-121, enero-junio. 2017. In the numerical simulation of the stochastic model it can be observed in Figure 2 that 6 different pathwise patterns occur with the mean trajectory included. Perturbations on b 0 in a range of 1% are used and Milstein numerical scheme [33] are applied to obtain the solutions of the stochastic system (3). All these trajectory solutions satisfy positivity and boundedness as the theoretical results predict. In order to give more support to our theoretical results, an extreme case with perturbations on b 0 in a range of 3% are used. As can be seen in Figure 3, the 10 pathwise trajectories with the mean trajectory are bounded and positive despite the use of large perturbations. It is important to remark that other perturbations in larger ranges have been taken and the positivity and boundedness are maintained.

|116
Ingeniería y Ciencia autores In Figure 4 and 5 we present 10 pathwise trajectories that show that the behavior of the infection tends to zero as t increases if the parameters satisfy β u < ν l . In addition, as it was expected in Figure 5 the susceptible proportion approximates to one as t increases. Thus, by computer simulation, we support our theoretical results for the conditions under the RSV die-out. Thus, health institutions can take measures to make the corresponding parameters change in order to obtain a feasible way for RSV prevention and control.

Conclusions
In this paper, we have considered a stochastic seasonal epidemic model in a population, where some stochasticity is introduced on the baseline transmission rate. Here, we investigated the existence and uniqueness of the solutions of the stochastic model and proved positivity and boundedness, which is of paramount importance for the study of the dynamics of population models. In addition, the positivity and boundedness of solutions is important to other nonlinear models that arise in sciences and engineering. Thus, a similar approach can be applied to other models from different areas.
Stochastic differential equations give another option to model viral dynamics and stochastic effects and introduce a more realistic way of modeling this type of disease. A numerical simulation of the seasonal stochastic models describing the transmission of the respiratory syncytial virus RSV in the region of Valencia using the Milstein numerical scheme is included in order to support our theoretical results. By computer simulation, we show how RSV die out when the a condition is satisfied. Thus, health institutions can take measures to make the corresponding parameter changes in order to obtain a feasible way for RSV prevention and control.