Mathematical modelling of respiratory syncytial virus (RSV): vaccination strategies and budget applications

Summary

We propose an age-structured mathematical model for respiratory syncytial virus in which children aged <1 year are especially considered. Real data on hospitalized children in the Spanish region of Valencia were used in order to determine some seasonal parameters of the model. Weekly predictions of the number of children aged <1 year that will be hospitalized in the following years in Valencia are presented using this model. Results are applied to estimate the regional cost of paediatric hospitalizations and to perform a cost-effectiveness analysis of possible vaccination strategies.

INTRODUCTION

Respiratory syncytial virus (RSV) is the most important respiratory virus in young children, and the major cause of hospitalizations, especially for bronchiolitis and pneumonia, in infants [1]. Its impact on health services is increasing as the incidence of hospitalization in children for bronchiolitis increases [2]. It is only recently that the impact of RSV on adults has been studied, as up to 18% of the pneumonia hospitalizations in patients aged >65 years are due to RSV [3].

In Spain, there are 15 000–20 000 visits to primary care due to RSV each year. In the Spanish region of Valencia, 1280 children aged <5 years are hospitalized each year as a consequence of RSV bronchiolitis [4]. The average hospitalization time is about 6 days. The cost of paediatric hospitalization for the Valencian Health Service is about €3·5 million per year.

RSV is the cause of annual seasonal epidemics with minor variations each year and its coincident occurrence with other widespread viral infections such as influenza or rotavirus, produces a high number of hospitalizations overstretching the health service. Moreover, it is easily transmitted and nosocomial infections are frequent [5].

Therefore, research on RSV and other viruses and the development of strategies to control epidemics are important from both the sanitary and economic point of view. Another problem is the study of vaccines to protect individuals at an early age, when the immune system is not completely developed. In this context it is important to note that the most seriously ill cases are not due to the RSV infection but are related to an anomalous immune response of the child [6].

Mathematical models have been shown to be a powerful tool to analyse the epidemiology of infectious illnesses, to understand their behaviour, to predict their social impact and to discover how external factors change the impact of disease. In the case of RSV, the building of a reliable model is a priority objective to predict the medical care requirements needed in the following seasons.

Mathematical models for RSV have been developed previously. For example, SIRS (susceptible–infectious–recovered–susceptible) and MSEIRS (maternally derived immunity–susceptible–latent–infectious–recovered–susceptible) mathematical models with four possible re-infections were studied and applied with data from The Gambia, Singapore, Florida and Finland by Weber et al. [7]. White et al. [8] presented a nested RSV model, stochastic simulations and fitting with data from several countries. Moreover, in an earlier paper, White et al. [9] considered the two types, A and B, of RSV and developed a SIRS model where re-infection by either type of RSV virus is possible, fitting the model with data from England & Wales and Finland.

Further, a PIV-vectored vaccine is already under development and clinical studies have been carried out recently [10, 11]. This vaccine might be available in the near future and, consequently, planning of vaccination strategies is urgently required. In this paper we consider a newborn vaccination strategy for Valencia in terms of the estimated cost of the vaccine, the average cost of hospitalization of RSV children that develop acute symptoms (in this case the average hospitalization time is ~6 days [4]) and parent work loss. To the best of our knowledge, vaccination strategies for RSV have not been studied and the imminence of the application of PIV-vectored vaccines demands such a study. A previous work on the cost-effectiveness of immunoprophylaxis with palivizumab has been recently reviewed [12].

Taking into account that RSV is an illness for which the most affected individuals (i.e. most hospitalizations) in the Spanish region of Valencia are children aged <1 year required the development of an adequate demographic model. We used a model with two age groups: G ₁ corresponding to children aged <1 year and G ₂ for the rest of the population. Birth and mortality rates as well as the transition rate from G ₁ to G ₂ were obtained from real statistical data.

The evolution of infected individuals is predicted by means of Weber et al.'s SIRS [7] generalized to take into account two age groups and numerically integrated by standard numerical methods [13].

In the following sections we discuss the age-structured SIRS model estimating some demographic parameters from data from the Instituto Valenciano de Estadística (Statistical Service of the Spanish region of Valencia) [14]. Data corresponding to hospitalizations of children aged <1 year are presented and a Nelder–Mead procedure for model fitting is developed. A cost-effectiveness analysis is performed for a vaccination scenario in the region of Valencia, Spain. Finally, conclusions are given.

METHODS

The Spanish region of Valencia [15] is located in eastern Mediterranean Spain, with an area of 23 255 km² and a population of 4 202 607 inhabitants (2001), composed of three provinces, Castellón (north), Alicante (south) and Valencia (middle).

Following [16, section 6·1, pp. 634–635], we introduced an age-structured SIRS model with two age groups: i=1 corresponds to children aged 0–1 year and i=2 corresponds to the rest of the population. This is justified because the disease is more acute in children aged <1 year, as confirmed by hospitalization data. For each age group we have three subpopulations according to the state of the individuals with respect to the disease:

• • Susceptibles, S _i(t) (i=1, 2): the proportion of those at risk of contracting the disease.

• • Infectives, I _i(t) (i=1, 2): the proportion of those infected and capable of transmitting the disease.

• • Recovered, R _i(t) (i=1, 2): the proportion of those recovered from the disease that are temporarily immune to re-infection.

Note that S _i(t), I _i(t) and R _i(t), respectively, correspond to the fraction of susceptible, infected and recovered of age group G _i (i=1, 2) and

The age-structured model is then defined by the following system of differential equations:

(1)

(2)

(3)

(4)

(5)

(6)

with the initial condition

(7)

(8)

where

• • μ is the birth rate;

• • β(t)=b ₀+b ₁ cos(2πt+ϕ), is the disease transmission rate, expressed as such in order to represent the seasonality of this disease. This type of function is also considered in other studies where seasonality is relevant [7, 17];

• • c ₁, is the transition rate from age groups G ₁ to G ₂;

• • d_i (i=1, 2), is the death rate for each age group;

• • γ_i, where 1/γ_i (i=1, 2) is the average time an individual of age group G _i remains immune against re-infection,

• • ν_i, where 1/ν_i (i=1, 2) is the average time to recover from the illness for an individual of age group G _i.

It should be noted that this model assumes some simplifying hypotheses:

1. (i) Primary and secondary infections have the same rate of recovery and infectivity but this is not the case in real clinical situations [5].

2. (ii) The mixing between the different age groups is homogeneous.

3. (iii) There is no maternal antibody protection.

It would be possible to generalize further the model in equations (1)–(6) to take into account that the strength of successive re-infections is progressively smaller. However, there is no reliable data for these clinical parameters for the different age groups. Heterogeneity could also be included by modulating the transmission rate, β(t), with a contact matrix C _ij (i, j=1, 2). But this would increase the number of parameters to fit and therefore, create an additional difficulty in obtaining reliable fitted parameters in an acceptable time. For these reasons, we decided to study the simpler model, taking into account that the cost-effectiveness of the vaccination programme could be slightly different in a more general model.

Maternal antibody protection has already been studied in the context of RSV models by Weber et al. [7], but there were no noticeable changes in respect of the simpler SIRS model without early antibody protection. On the other hand, we performed CHAID (χ² automatic interaction detection) analysis with data of hospitalized children in the Spanish region of Valencia, obtaining an average duration of 23 days for the maternal antibody protection period, and this can be considered as a perturbation of the simpler SIRS model.

Estimating parameters

Numerical values can be assigned to most of the above parameters. For instance, γ_i, and ν_i (i=1, 2) can be obtained from [1, 7], the average time to recover from RSV is 10 days and the average time to lose immunity is 200 days, both independent of the age group. Then,

(9)

(10)

If we denote by N _i the fraction of persons in age group G _i (i=1, 2) (in respect of the total population) the evolution equations for the fraction of individuals in each group is given by:

(11)

(12)

(13)

Using demographic data of the Spanish region of Valencia for the period 2001–2004 [14], where hospitalization data are available, we obtain:

• • the mean population, 4 252 386 inhabitants;

• • the mean fraction of both subpopulations, N ₁=0·01038 and N ₂=0·98962;

• • the mean birth rate, μ=0·0107497;

• • the mean death rate for each age group, d ₁=0·0006796 and d ₂=0·00912862.

Under the assumption of constant population for each age group (N ₁^′=N ₂^′=0) and equation (11), we find a simple relation for the transition rate from age group G ₁ to age group G ₂:

(14)

By adding equations (11)–(13) we also find the following relation:

(15)

Finally, we choose a new d ₂ to satisfy equation (15), i.e. d ₂*=0·0108554, which is different from the real value of d ₂=0·00912862 because the population is not really constant in the 4-year interval considered. Nevertheless, the difference is sufficiently small (~7200 dead persons) to be ignored in our study.

Hence, as can be seen, only parameters b ₀, b ₁ and ϕ remain unknown, i.e. parameters dependent upon the seasonality.

Source RSV data

From CMBD (basic minimum dataset) database of the Spanish region of Valencia we obtained hospitalization data on RSV-related illness (bronchiolitis, infection, pneumonia, etc.) of children aged <1 year during the period January 2001 to December 2004. After processing the data we obtained the weekly number of hospitalizations of children in age group G ₁ caused by a RSV-related infection as shown in Figure 1.

Fig. 1. Weekly hospitalizations for respiratory syncytial virus (RSV) infection from January 2001 to December 2004 of children aged <1 year in the region of Valencia.

It should be noted that data refer to hospitalizations, not to infected people. This means that only a proportion of infected are hospitalized and this proportion or scale, designated s, should be estimated in the model data fitting.

RESULTS

Model fitting

In order to fit the data for hospitalized children aged <1 year we numerically integrated the system of differential equations in equations (1)–(5) for different sets of unknown parameters: b ₀, b ₁, ϕ and the scale parameter s of the proportion of infected in G ₁ that are hospitalized.

Additionally, the initial condition data are also unknown. We have considered throughout all the simulations that 1% of children in age group G ₁ are infected at 1 January 1999 and the rest of individuals are susceptible.

The procedure for fitting the model is as follows: we solve numerically the system of differential equations in equation (1)–(5) with the aforementioned initial values. After a period of 2 years of transient evolution the solution becomes periodic. This numerical solution is optimized by mean of the Nelder–Mead algorithm [18, 19], that does not need the computation of any derivative or gradient, which is impossible to know in this case.

The values of b ₀, b ₁, ϕ and s that minimize the function for the region of Valencia are:

(16)

(17)

(18)

(19)

For the different health districts where the fitting procedure was performed we obtained very similar values. This supports the reliability of the fitting procedure because the parameters b ₀, b ₁ and ϕ are known to depend only on the climatological conditions of the countries where the RSV epidemic has been studied [8, 9].

The parameter s is especially interesting because it provides information about the proportion of infected children that have to be hospitalized as a consequence of the infection. According to the data in equation (19), we found that 2·19% of infected children are hospitalized by a disease derived from RSV infection.

The graphical representations of this model fitting can be seen in Figure 2.

Fig. 2. Model fitting from January 2001 to December 2004. The points are data of weekly hospitalizations of children aged <1 year and the continuous line the scaled solution of the system of differential equations for infected children aged <1 year, s*I ₁(t). The dashed line corresponds to the fitting with γ₁=γ₂=1·59.

The average time of recovery from RSV illness (ν₁=ν₂=365/10) is a well-known value [1, 7]. It can also be obtained and compared with the statistical analysis of hospitalization data.

On the other hand, a larger uncertainty is found in the average time of immunity following RSV infection, 1/γ. Weber et al. [7] reported the value γ=1·8 corresponding to an average of 200 days of immunity. This value was contrasted by a sensibility analysis. Similar values have been proposed by White et al. [8]. In order to verify that the fitting of RSV hospitalization data to a SIRS model requires and immunity period compatible with the values suggested in the literature we performed a more detailed sensibility analysis. We considered that the value of γ is in the range 1·0<γ<2·5. A Nelder–Mead optimization algorithm was applied to obtain the best values of γ, the seasonal parameters b ₀, b ₁, ϕ and the scale factor s. The best fit is obtained with the values:

(20)

(21)

(22)

(23)

(24)

These values are very similar to those obtained in equations (16)–(19) with the assumption γ=1·8 and confirm that a SIRS model approach to the spread of RSV disease implies an immunity time after infection of ~200 days. In Figure 2 we also compared this with a second fitting; however, the difference was minimal.

Vaccination strategies

We next developed a cost-effectiveness analysis for a newborn vaccination strategy for RSV. The effect of the vaccine is to remove a fraction p of newborns from the S ₁ state and place them on a new vaccinated V state until they reach the age of the next age group. The box diagram for the model is plotted in Figure 3.

Fig. 3. Box diagram of a respiratory syncytial virus model with two age groups and newborn vaccination.

The reason that V state individuals evolve to the S ₂ state instead of S ₁ is the large immunization period, which is at least equal to the immunization of children recovered from the disease (γ=365/200 days). Moreover, the children should receive additional booster doses during the first year of life and, consequently, their immunization is prolonged.

Figure 4 plots the evolution of the total number of infected children aged <1 year in the region of Valencia after implementation of the vaccination programme immediately after the beginning of the fifth year. A total of 85% is vaccinated. It can be seen that the number of infected children decays to <1000 just 2 years after commencement of vaccination.

Fig. 4. Prediction of the model for the total number of infected children in the region of Valencia. Vaccination of 85% of newborn children begins after the fifth year.

The total cost of RSV healthcare is calculated as follows:

• • Hospitalization cost: Taking into account an average of 6·28 hospitalization days for every acutely infected child [4] and €500 per day and child hospitalized.

• • Vaccination cost: We calculated the number of children vaccinated during a year and an estimated cost of €100 per dose. Three doses are planned during the first year of life.

• • Parent work loss: We take into account that for every hospitalized child parents lose d=6·28 days of work (average hospitalization of RSV-infected children). In the case of infected children we assume that they develop milder symptoms but, nevertheless, parents lose d days of work. This is a social parameter difficult to determine. We simulated cost-effectiveness for d=2, 3, 4. The mean labour cost in Spain is €75·21 per day [20].

The results for the costs of hospitalization, vaccination and parent work loss are plotted in Figures 5 and 6.

Fig. 5. Vaccination cost (—) and hospitalization cost (……) vs. the fraction of vaccinated newborn children.

Fig. 6. Costs for parent work loss. Taking into account d=6·28 days of work loss on average for hospitalized children and d=4 (—), d=3 (- - -) and d=2 (……) days of work loss for infected children that do not need hospitalization.

The global result for the cost-effectiveness analysis is plotted in Figure 7 for d=2, 3 or 4 days of parent work loss in the case of children that do not develop sufficiently acute symptoms to become hospitalized. A reduction of €2 million of total cost is predicted for an estimate of 4 days of parent work loss on average for infected children. Moreover, we have not taken into account the long-term effects of RSV infections. In particular, there is an accord among paediatricians regarding a connection between RSV at early ages and asthma episodes of children and adolescents. This has been confirmed by recent studies in mice [21]. Therefore, even assuming that parents only lose four working days in caring for children who develop mild symptoms of RSV, a positive balance for the implementation of the vaccine is obtained.

Fig. 7. Total cost (hospitalization+parent work loss+vaccination) related to respiratory syncytial virus infections as a function of the fraction of vaccinated newborn children. Days of parent work loss: d=4 (—), d=3 (- - -) and d=2 (……).

DISCUSSION

In this paper we have presented an age-structured model for RSV that fits well with the data of hospitalized children aged <1 year in the region of Valencia. A Nelder–Mead method [18] has been implemented in order to estimate the parameters which determine the seasonality and also the fraction of infected children that become hospitalized as a consequence of the RSV infection.

We used data of children aged <1 year hospitalized as a consequence of RSV infection in the period 2001–2004.

Independently of the health district considered, we have obtained very similar parameters for the variation of the infection rate with time in the SIR model: β(t)=b ₀+b ₁ cos(2πt+ϕ). This is unsurprising as it is expected that these parameters depend only on the climatological conditions of the country [8, 9]. The climate in Valencia is Mediterranean with an average temperature around 18°C and minimum temperatures in January (the month in which the highest number of child hospitalizations due to RSV is observed). Moreover, the absence of mountain ranges, except for the Iberian Mountains in the northwest where population density is small, implies no significant climatic changes between the coast and the interior of the region. The study of these correlations is an interesting problem in itself and must be treated by means of statistical analysis. Sociodemographic factors, e.g. the opening and closing of schools, could be a secondary source for the outbreak of RSV. However, this is apparently not as important as climate because in Valencia schools open in September and the seasonal outbreak normally occurs at the end of November.

We have developed a model for newborn vaccinations with two booster doses during the first year of life. Assuming a cost of €300 per vaccination a positive balance for the implementation of the proposed vaccination strategy is obtained even with the cautious assumption that parents only lose 4 working days in the case of an infected child that does not become hospitalized. We must take into account that these results are obtained with a simple extension of a recently proposed SIRS model [7] and should be considered as a first-order approximation. However, we are also aware that a mean-field network approach replicates the results discussed in our paper even with a continuous Forster–McKendrick population model [16], which is more realistic than the two age groups proposed in our paper.

Another possible alternative is the MSEIRS model with re-infections [7]; however, taking into account that the vaccine is expected to prevent infections during the first year of life similar results to the model proposed in our paper are expected.

On the other hand, although the cosine function to simulate the seasonal strength of the infection rate, β(t), is the most widely used, other alternatives utilized in relation to climatological factors are probably more realistic. Nevertheless, there is no general agreement on which factors are relevant in unleashing the most acute episodes of the epidemic during the winter season: humidity, temperature or even ultraviolet B radiation has been linked with the propagation of RSV virus [22–24]. Further research on these factors in the region of Valencia could also provide better information on the epidemic propagation and the outcome of vaccination programmes. We are currently performing this statistical correlation analysis.

A more general model including heterogeneity of contacts in individuals, variable strength of re-infection, a more realistic model for β(t) based upon climatological factors and maternal antibody protection is currently under study and will be published elsewhere.