Parameter Estimation, Sensitivity Analysis and Optimal Control of a Periodic Epidemic Model with Application to HRSV in Florida

A state wide Human Respiratory Syncytial Virus (HRSV) surveillance system was implemented in Florida in 1999 to support clinical decision-making for prophylaxis of premature infants. The research presented in this paper addresses the problem of fitting real data collected by the Florida HRSV surveillance system by using a periodic SEIRS mathematical model. A sensitivity and cost-effectiveness analysis of the model is done and an optimal control problem is formulated and solved with treatment as the control variable.


Introduction
Human respiratory syncytial virus (HRSV) is a virus that causes respiratory tract infections. It is a major cause of lower respiratory tract infections and hospital visits during infancy and childhood. A prophylactic medication, palivizumab, can be employed to prevent HRSV in preterm infants (under 35 weeks gestation), infants with certain congenital heart defects or bronchopulmonary dysplasia, and infants with congenital malformations of the airway. Treatment is limited to supportive care, including oxygen therapy. In temperate climates, there is an annual epidemic during the winter months. In tropical climates, infection is most common during the rainy season. In the United States, 60% of infants are infected during their first HRSV season, and nearly all children will have been infected with the virus by two to three years of age [9]. Of those infected with HRSV, 2 to 3% will develop bronchiolitis, necessitating hospitalization [10]. Natural infection with HRSV induces protective immunity, which wanes over time, possibly more so than other respiratory viral infections, and thus people can be infected multiple times. Sometimes an infant can become symptomatically infected more than once, even within a single HRSV season. Severe HRSV infections have increasingly been found among elderly patients. Young adults can be reinfected every five to seven years, with symptoms looking like a sinus infection or a cold. The Florida Department of Health provides an integrated and reliable HRSV system, with data from hospitals and laboratories [8].
Mathematical models can project how infectious diseases progress, to show the likely outcome of an epidemic, and help inform public health interventions. In epidemiology, compartmental models serve as the base mathematical framework for understanding the complex dynamics of these systems. Such compartments, in the simplest case, stratify the population into two health states: susceptible to the infection of the pathogen, often denoted by S, and infected by the pathogen, often denoted by the symbol I. The way that these compartments

Review of some recent HRSV mathematical models
We focus on compartmental models that divide the population into mutually exclusive distinct groups (of susceptible, or infected, or immune individuals, or . . . ) and we use deterministic continuous transitions between those groups, also known as states. Due to the seasonality of HRSV, the models that best fit real data are periodic and, therefore, non-autonomous. In [30], two models are proposed where the transmission is periodic: (i) a simple model with only three compartments, known as SIRS (Section 2.1); (ii) and a more complex model with seventeen compartments, named MSEIRS4. However, it is shown that the simpler model is closer to real data [30]. Zang et al. [31] use a non-autonomous SEIR model where, beyond the periodicity in the transmission rate, the annual recruitment rate is also periodic. This assumption is due to opening and closing of schools [31]. Here we adopt such ideas and consider a simple non-autonomous and periodic SEIRS model (Section 2.2).

SIRS model
In the SIRS model, one considers that the population consists of susceptible (S), infected and infectious (I), and recovered (R) individuals. A characteristic feature of HRSV is that immunity after infection is temporary, so that the recovered individuals become susceptible again [30]. The model depends on several parameters: parameter µ, which denotes the birth rate, assumed equal to the mortality rate; γ, which is the rate of immunity loss; and ν, which is the rate of loss of infectiousness. Moreover, the influence of the seasonality on the transmission parameter β is modeled by the cosine function. Precisely, and using the linear mass action law [30], the SIRS model for HRSV is given by the following system of ordinary differential equations [30]: where β(t) = b 0 (1 + b 1 cos(2πt + Φ)). Note that b 0 is the average of the transmission parameter β, b 1 is the amplitude of the seasonal fluctuation in the transmission parameter β, while Φ is an angle that will be chosen later in agreement with the real data for I(t).

SEIRS model
In order to incorporate some important features of HRSV into the model while keeping it simple, we extend the model (1) of Section 2.1 by using the ideas in [13]: firstly, we include a latency period by introducing a group E of individuals who have been infected but are not yet infectious, these individuals becoming infectious at a rate ε and the latency period being equal to the time between infection and the first symptoms; secondly, and similar to [31], we consider that the annual recruitment rate is seasonal due to schools opening/closing periods. The equations of the SEIRS model are then given by where λ(t) = µ(1 + c 1 cos(2πt + Φ)) is the recruitment rate, including newborns, immigrants, etc. Parameter c 1 is the amplitude of the seasonal fluctuation in the recruitment parameter λ and Φ an angle to be conveniently chosen.

Main results
We begin by investigating how realistic the models discussed in Section 2 are with respect to HRSV and real data from Florida [8]. For that we need a proper estimation of the parameter values.

Parameter estimation
During model fitting, the values of the parameters µ (birth rate), ε (infectious rate), ν (rate of infectiousness loss) and γ (rate of immunity loss) were held constant. See Table 1, where we use n to denote the average number of monthly HRSV cases reported. The values of ε, ν and γ were taken from [30]. For the birth rate, µ, we used the value given in [7] for the state of Florida. In agreement with Section 2.1, we have set the birth rate equal to the mortality rate, so that one obtains a constant population during the period under study. This simplification is justified when, as in this case, the annual infection rate is much bigger than the growth of population. To simplify analysis, we set the total population equal to one and use fractions for the values of state variables. Values of four parameters were determined by fitting the model: (i) the mean of the transmission parameter, b 0 ; (ii) its relative seasonal amplitude, b 1 ; (iii) the angle Φ; and (iv) a scaling factor s, which scales the number of infectious individuals to the empirical case reported in our unit scaled model. The angle Φ was normalized in the following way: for a value of zero, a maximum of the cosine function used in β(t) (and also used in λ(t)) coincides with the first maximum of the empirical cases. The models were fitted to data on the reported number of positive tests of HRSV disease, per month, in the State of Florida, excluding North region, between September 2011 and July 2014 (35 months). The data was obtained from the Florida Department of Health [8]. Our results are given in Figure 1.  Following [28], our fitting approach consisted to minimize the l 2 norm of the difference between the real data and predictive cases of HRSV infection given by models (1) and (2). Let e represent the relative error, computed by  Table 1. where I empiric is the real data obtained from [8] and I model the one predicted by the SIRS or SEIRS models. The results shown in Figure 1 correspond to a relative error of 0.066% and 0.060% of infants per year with respect to the total child population of Florida in 2014, excluding North region, respectively for SIRS or SEIRS models. We note that the values of µ, ν, γ and ε are annual rates. The difference between the two models, in terms of absolute errors, clearly justifies SEIRS' superiority.

Sensitivity analysis
One of the most important thresholds while studying infectious disease models is the basic reproduction number [29]. The basic reproduction number of the SIRS model was computed in [30]. Assuming that the transmission and the recruitment parameters are constant, the basic reproduction number for the SEIRS model is given by (see [16]) Now we do a sensitivity analysis for the basic reproduction number (3). Such analysis tells us how important each parameter is to disease transmission. This information is crucial not only for experimental design, but also to data assimilation and reduction of complex nonlinear models [20]. Sensitivity analysis is commonly used to determine the robustness of model predictions to parameter values, since there are usually errors in data collection and presumed parameter values. It is used to discover parameters that have a high impact on R 0 and should be targeted by intervention strategies. More precisely, sensitivity indices's allows to measure the relative change in a variable when parameter changes. For that we use the normalized forward sensitivity index of a variable, with respect to a given parameter, which is defined as the ratio of the relative change in the variable to the relative change in the parameter. If such variable is differentiable with respect to the parameter, then the sensitivity index is defined using partial derivatives, as follows (see [4,22]).

Definition 1
The normalized forward sensitivity index of R 0 , which is differentiable with respect to a given parameter p, is defined by One can easily compute an analytical expression for the sensitivity of R 0 , using the explicit formula (3), to each parameter that it includes. The values of the sensitivity indices for the parameters values of Table 1 are presented in Table 2. Note that the sensitivity index may depend on several parameters of the system, but also can be constant, independent of any parameter. For example, Υ R0 β = +1, meaning that increasing (decreasing) β by a given percentage increases (decreases) always R 0 by that same percentage. The estimation of a sensitive parameter Table 2. Sensitivity of R 0 evaluated for the parameter values given in Table 1, according with (3) and Definition 1.

Parameter
Sensitivity should be carefully done, since a small perturbation in such parameter leads to relevant quantitative changes. On the other hand, the estimation of a parameter with a small value for the sensitivity index does not require as much attention to estimate, because a small perturbation in that parameter leads to small changes [14]. According with Table 2, we should pay special attention to the estimation of parameters ν, b 0 (mean value of β(t)) and ε. In contrast, the estimation of the rate of birth, µ, (or rate of mortality) does not require as much attention because of its low value of the sensitivity index. This is well illustrated in Figure 2, where we can see in Figure 2a the graphics of the number of infectives with and without an increment of 10% for the parameter ν (most sensitive parameter), and in Figure 2b for parameter µ (less sensitive parameter).   Table 1 and with an increase of 10% of a specific parameter: ν in (a) and µ in (b).

Optimal control
To investigate some optimal control measures, we choose the SEIRS model, which provides a better fitting of available real data, as shown in Section 3.1. The evolution on the number of susceptible, exposed, infective and recovered, depend on some factors that can be controlled. In case of HRSV, treatment is the most realistic. For this 142 PARAMETER ESTIMATION, SENSITIVITY ANALYSIS AND OPTIMAL CONTROL OF A PERIODIC EPIDEMIC MODEL reason, we consider the following optimal control problem:

t)S(t)I(t) + γR(t), E(t) = β(t)S(t)I(t) − µE(t) − εE(t),
subject to given initial conditions S(0), E(0), I(0), R(0) 0, where Ì denotes treatment (Ì is the control variable). Note that in the absence of control measures, that is, when Ì(t) ≡ 0, system (4) reduces to (2). We consider the set of admissible control functions Our purpose is to minimize the number of infectious individuals and the cost required to control the disease by treating the infected. Precisely, the optimal control problem consists in with 0 < κ 1 , κ 2 < ∞. To solve the problem, we apply the celebrated Pontryagin maximum principle [19]: the Hamiltonian is given by the maximality condition ensures that the extremal control is given by while the adjoint system asserts that the co-state variables p i (t), 1 ≤ i ≤ 4, satisfy relationṡ and transversality conditions Therefore, in order to solve the optimal control problem (6), we solve the following boundary value problem: with given initial conditions and transversality conditions (9), where Ì(t) is given by (7). This is done numerically in Section 3.4.

Numerical results and cost-effectiveness analysis for the optimal control model
Two algorithms were implemented to obtain and confirm the numerical results. One approach solves the optimal control problem numerically using a fourth order Runge-Kutta iterative method. First we solve the system (4) with initial conditions for the state variables and a guess for the control over the time interval [0, t f ], by the forward Runge-Kutta fourth order procedure, and obtain the values of the state variables S, E, I and R. Using those values, then we solve the system (8) with the transversality conditions (9), by the backward fourth order Runge-Kutta procedure, and obtain the values of the co-state variables. The control is updated by a convex combination of the previous control and the value from (7). The iteration is stopped when the values of the unknowns at the earlier iteration are very close to the ones at the current iteration. The results coincide with the ones obtained by a similar iterative method that, in each iteration, solves the boundary value problem (10) using the MATLAB bvp4c routine.
In what follows, we consider that, for simplicity, Ì max = 1 and the other parameters are fixed according to Table 1, with exception to the angle Φ that is assumed to be π/2. Such value allows that the transmission parameter initial value be the average, β(0) = b 0 , and the recruitment rate initial value also be the average, λ(0) = µ. The initial conditions, given by Table 3, are obtained as the nontrivial equilibria for the system (4) with no control, corresponding to the population state prior the introduction of the treatment. We also assume that t f = 5 because the World Health Organization goals for most diseases are usually fixed for 5 years periods. Table 3. Initial conditions for the optimal control problem with parameters given by Table 1, excepting angle Φ that is assumed to be π/2. The values correspond to the endemic equilibrium of (4) before the introduction of treatment.
The solution for the optimal control problem is illustrated in Figures 3, 4a and 4b. The periodic nature of the disease influences the evolution of the four state variables. We can also see that the control is a continuous function that shows some nonregularity at the end of the time interval  levels. Figure 4c shows the efficacy function [24] that is defined by where I * (t) is the optimal solution associated with the optimal control and I(0) is the corresponding initial condition. This function measures the proportional variation in the number of infectious individuals after the intervention with control Ì * , by comparing the number of infected individuals at time t with the initial value I(0) for which there is no control implemented. We observe that F (t) oscillates between −1.36 (lower bound) and +0.68 (upper bound), and exhibits the inverse tendency of I(t).
Obviously, the results depend on the objective functional J given in (6). Namely, they depend on the weight constants associated with the amount of infectious individuals k 1 and with the cost of the control k 2 . According with Figure 5, results do not change qualitatively by varying constants k i , i = 1, 2. Nevertheless, the magnitude of the efficacy changes slightly when k 1 and k 2 vary independently.
Some summary measures are introduced to evaluate the cost and the effectiveness of the proposed control measure for the intervention period.
The total cases averted by the intervention during the time period t f [24] is given by where I * (t) is the optimal solution associated with the optimal control Ì * and I(0) is the corresponding initial condition. Note that this initial condition is obtained as the equilibrium proportion I of system (4)   intervention, which is independent on time, so t f I(0) = t f 0 I dt represents the total infectious cases over a given period of t f years.
We define effectiveness as the proportion of cases averted on the total cases possible under no intervention [24]: The total cost associated with the intervention [24] is where C corresponds to the per person unit cost of the detection and treatment of infectious individuals. Following [18,24], the average cost-effectiveness ratio is defined by Table 4 summarizes the particular case we have analyzed. These results evidence the limited effectiveness of the control variable to diminish the HRSV infectious individuals. Efficacy function F(t)

Conclusion
Human Respiratory Syncytial Virus (HRSV) is the most common cause of lower respiratory tract infection in infants and children worldwide. In addition, HRSV causes serious disease in elderly and immune compromised individuals. In this work, we discussed a mathematical compartmental model for HRSV. Estimation of parameters was done for real data of Florida from September 2011 to July 2014, minimizing the l 2 norm. The results show that the proposed model fits well the reality under study. Moreover, a sensitivity analysis was carried out for the basic reproduction number, in the case when the transmission parameter is taken to be the average value of β(t) in the period under study, proving that the most important parameters to have into account are the rate of loss of infectiousness ν and the average of transmission parameter b 0 . Our results from optimal control show that treatment has a limited effect on HRSV infected individuals. This reinforces the importance of developing a licensed vaccine for HRSV, which is a subject under strong current development [17].