Modeling the Dynamic of Multi-Wave Diseases: The Model of Hand, Foot and Mouth Disease

1. Introduction

Hand, foot and mouth disease (HFMD) is often a multiphase disease involving a succession of different viruses and waves. The co-circulation of different serotypes and alternation between enterovirus A71 (EV-A71) and coxsackievirus A (CV-A ) was commonly observed during HFMD epidemics [1,2,3,4]. However, co-infection is not a primary cause of severe forms [5,6]. Epidemics could be due to both the accumulation of susceptible individuals in the community, patients are mostly young children under five and introduction of new genotypes or strains [7]. Within the West Pacific – South East Asia region, outbreaks in Taiwan in 1998 and 2000 were caused by EV-71 C2 and B4 strains, respectively [8]. Sentinel surveillance in Sarawak, Malaysia, demonstrated that the emergence of the novel subgenotype C1 of EV-71 was the cause of the 2003 outbreak [9]. Although very often mild, HFMD can result in severe complications such as encephalitis, aseptic meningitis, pulmonary edema, myocarditis, and death [10] and can be devastating like in Cambodia in 2012 with a death rate of 88% [11]. This large panel of symptoms and variation over time depend for a good part on the succession of viruses driving the epidemics. Since this early period, incidence rate of HFMD in the Asia-Pacific region has been increasing, and China is the main endemic area [12].

A relevant issue in HFMD epidemiology could be to understand and distinguish the different waves during multipeak epidemic. Waves may be associated to different virus, with different transmissibility and virulence characteristics. Understanding the process of emergence of a disease can be related to the ability to identify type samples from outliers. Dynamic models are therefore required for analyzing epidemics. However, if current models can analyze individual waves, there is still a need to model and compute a multi-wave epidemic. The 2011-2012 HFMD epidemics was the largest experienced in Vietnam and the very first one to occur in North Vietnam. Hải Phòng displayed the highest prevalence in North Vietnam providing thus a large cohort of more than 9000 patients to analyze (Supplementary Table S1). Being the first HFMD epidemic in North Vietnam, this outbreak provided an opportunity to study the dynamic without interference from previous outbreaks. We thus report here a mathematical approach for discrimination of individuals within an outbreak time frame and a mathematical model for characterizing multi-wave outbreaks displaying differing wave-associated with virulence or transmissibility.

2. Material and Methods

2.1. Epidemiological Information and Specimen Collection

All HFMD cases in Hải Phòng province were reported to the National Institute of Hygiene and Epidemiology (NIHE) through the national communicable disease surveillance system since 2011.

2.2. Statistical and Mathematical Analysis

Mean comparison was implemented by a Student’s T-test. A Chi-square test was used for proportion comparison of Hải Phòng city population and a one-way ANOVA test was used for the variance analysis. The EMD distance was calculated using the EMD() function in Matlab. GIS analysis was performed using QGIS.

2.3. Bias and Ethics

Training session HFMD cases definition and reporting were organized for the staff of the routine surveillance system to enhance quality and consistency of case report. This work was conducted following the requirements of the Vietnamese Ministry of Health and under the Law of Communicable Diseases Prevention and Control passed in 2007.

3. Results

3.1. HFMD Burden during the 2011-2012 Epidemic

A large HFMD epidemic hit North Vietnam in 2011-2012 (65,949 cases, source: General Department of Preventive Medicine, Ministry of Health). It was the first outbreak recorded in this part of the country. The region of Hải Phòng was the hardest hit among the 28 North Vietnam provinces (Figure 1a). The epidemic in Hải Phòng province was slightly delayed in 2011 when compared to the rest of North Vietnam. The index case was a 6-years old girl admitted in an urban district hospital in week 16 (April 17, 2011). In total, 9621 cases were collected that could be of three waves during the time of the epidemics (Figure 1b & Supplementary Table 1). Several changes happen during the period of study that must be clarified. A first guideline concerned surveillance, prevention and control of HFMD was published on the February 24, 2012. The second guideline about diagnosis and treatment was issued on March 30, 2012. The publication these two specific guidelines by the Ministry of Health was necessary to improve the management of patients with severe signs of the disease. The influx of patients of all types to hospital led to overcrowding of the healthcare system during the first phase. The health system was reorganized during the second phase of the epidemic to deal only with the most severe cases. The result was an inversion of the curves between severe cases and those of lesser severity (Figure 1b). The etiological agent is a second element that needs to be taken into account in the overall picture of this epidemic in the Hải Phòng region. The first two phases of the epidemic were mostly associated with enteroviruses, mainly EV-71 according to molecular analysis performed on 257 throat swabs collected at the main pediatric hospital in Hải Phòng city. Wave 3 was associated with the co-circulation of CV-A6 and CV-A16.

3.2. Spatial Heterogeneity of Hải Phòng District

The spatial organization of the Hải Phòng province could have played a major role in the dynamic of expansion and evolution of the disease. HFMD transfer ways depend both on natural and anthropogenic environment parameters. The main factors are rivers, roads and settlements. The main West-East orientation of the rivers limits the North-South circulation whereas West-Est transit between coast and mainland is important (Figure 2a). Main road networks are located in distance from the main rivers, but follow the same main West-East orientation. North-South connections are less common and correspond to few secondary roads and rivers. The density of secondary roads increases strongly in the Hải Phòng city suburb area. Settlements are generally isolated from the main rivers, mostly because of flood risk. Rural areas are connected to settlements by province roads yielding a good locally connected settlement network, but they remain fragmented as pockets isolated by the main rivers from neighboring territories.‏ Hải Phòng city and its suburb area sprawl West-East without crossing the river to the north.‏ Settlement clusters south of Hải Phòng are regularly spread whereas on the northern part of the province high concentrations are found along the main country road while villages are scattered over the northern territory (Supplementary Figure 1). Correspondence of communes, ID numbers, district names and commune names are given in Supplementary Table S2.

The early cases appeared in the northern and central part of the region i.e., Hải Phòng city, (Supplementary Figure 2). First cases emerged in the district of Ngô Quyền before spreading south and west and north west. The southern part of the province occurrence of subsequent cases was strongly delayed. Only one early location was found in the southern part, corresponding to the largest urban area, while two locations are identified in southeastern coastal region. Early cases were also recorded on some commune of the main island. The main axis of virus expansion in the North is primarily West-East, spreading over the northern part of the province along rivers and the main country roads. The second wave started in the more southern district of Dương Kinh and then spread south and west before going back to the east. Finally, the third wave started in the central district of Kiến An and followed a similar route as wave 2.

3.3. Local Heterogeneity of the Epidemic

A film summarizes the number of cases over time in each commune (Supplementary video, https://filesender.renater.fr/?s=download&token=d75a1140-63fc-4df8-bb9a-ab20f19ee067). Foci of infection emerged sporadically. The frequency of variations in each commune increases during the epidemic waves. Some nearby communes appear to behave in a coordinated fashion or just spread out over time, suggesting possible transmission routes for the virus in the region. The EMD (Earth Mother Distance) method was used to compare case distribution between communes. The number of cases was reported by month for each commune to get enough data in low density commune. A distance matrix of all communes was used to build a tree and a classification based on 21 clusters that were then mapped on the communes (Figure 2b). Virus circulation followed some communication routes. The red communes (cluster 1 and 2) in the center of the region are along west part of QLB5 main road and a secondary route from main crossroad in the Hồ Nam commune of Lê Chân district. They correspond to the east to west transmission of the virus from the sites of emergence of waves 2 and 3. Cluster 2 was also found in the south-west part of the region for communes located along QL10 road. District level confirmed the heterogeneity of the epidemic in the region (Supplementary Figure 3). Most districts show three epidemic waves. The intensity of each wave varies from district to district. The most populous and popular urban districts show a distribution similar to that observed across the region as a whole. There is also a difference in patient management and adherence to guidelines published in February and March 2012 in one district in the south of the region and one urban district in the east. Local heterogeneity and changes of guidelines question the possibility to model the epidemic in three waves with a unique model.

3.4. Model Development

The in-depth characterization of the epidemic showed that within the one-year period of the study, three different waves occurred, each separated by a minimum, i.e. a drop in the number of cases followed by an increase. This in depth-analysis also showed that the diversity of symptoms is not a criterion for characterizing a wave as it is dependent on parents and physician behavior who apply a precautionary attitude and exaggerate the symptoms to ensure an immediate care. To find objective parameters each step was thus analyzed separately using the same Bernoulli model. In a first step, we assumed that each wave is relatively well separated. It means that when a wave is closing when the next wave is starting in a sense to be specified later. The function was therefore aiming at describing HFMD cases over time. The total number of HFMD cases for one wave is N. Over a defined period called h, the number of new cases depends of the number of susceptible patients in the population Xs(t) and two parameters related to the epidemic: the probability K(t) to be exposed to the virus and K0 the probability to present clinical symptoms. These two probabilities reflect at epidemiological level two the virus associated parameters which are virulence and pathogenicity, respectively.

The cumulative number Xm(t+h) of patients, is governed by the equation for one wave:

$$X_m(t+h)=X_m\left(t\right)+K_{0}K(t)X{}_S\left(t\right)*h,t\ge t_0\left(1\right)$$

with:

$$X_s\left(t\right)=N-X_m\left(t\right)\mathrm{and}K\left(t\right)=\frac{X_m\left(t\right)}{N}$$

The dot K(t) Xs(t) is the expected number of exposed patients. The constant K0 is called interaction factor. This is the probability of an exposed patient to become sick.

By passing to the limit as h tends to 0, we obtain a differential equation of Bernoulli (2)

$$\frac{X_m^{\text{'}}\left(t\right)}{X_m^2\left(t\right)}-K_0\frac{1}{X_m\left(t\right)}+\frac{K_0}{N}=0,\forall t\ge t_0\mathrm{and}{X}_m\left(t\right)\ne 0\left(2\right)$$

with the two limits conditions:

The general solution of equation (2) is the well-known logistic function defined by three parameters (K0, K1, K2):

$${\mathit{X}}_{\mathit{m}}\left(\mathit{t}\right)=\frac{{\mathit{K}}_2}{(1+\mathbf{e}\mathbf{x}\mathbf{p}(-{\mathit{K}}_1-{\mathit{K}}_0(\mathit{t}-{\mathit{t}}_0)},\forall \mathit{t}\ge {\mathit{t}}_0$$

(3)

With the help of the vector Y of observed values ${\mathit{Y}\left({\mathit{t}}_{\mathit{i}}\right)}_{0\le \mathit{i}\le \mathit{k}}$, the constant K0, K1, K2 are computed with a jacquard algorithm by minizing

$$Ec\left(K_0,K_1,K_2\right)={\displaystyle \sum _{i=0}^{\tau }{\left(y\left(t_i\right)-X_m\left(t_i\right)\right)}^2}$$

(4)

Note that τ_κ is the experimental end of the wave. This corresponds to the time of the last observed case of a given wave. The time t0 is the date of the first case observed for a given wave. Nevertheless, equation (3) can be interpreted with two limits conditions:

$$\underset{t\to +\mathrm{\infty }}{\mathrm{Lim}}Xm\left(t\right)=NandXm\left(t0\right)=N_0$$

N and N0 are theoretical values. They are respectively the total number of cases at the theoretical end of the waves $t=+\infty$ and the number of cases at the beginning of the wave. With equation (3), we have N=K2 and $N_0=\frac{K_2}{(1+\mathrm{e}\mathrm{x}\mathrm{p}(-K_1)}$).

Owing to the limit conditions and N>N0, it is possible to demonstrate that the interaction factor K0 is positive and less than 1. Therefore, this factor can be seen as a probability.

3.5. Wave-Specific Model

We assumed that the time of the beginning of each wave (t0) is known. Therefore, the beginning of the wave 2: ${\mathrm{t}}_{\mathrm{0,2}}$ occurred after the date of the last observed case of wave 1 ${\mathsf{\tau }}_1$. The same condition is assumed for wave 2 and wave 3. This means that: ${\mathrm{t}}_{\mathrm{0,2}}>{\mathsf{\tau }}_1\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{t}}_{\mathrm{0,3}}>{\mathsf{\tau }}_2$(Figure 3a). These dates are provided by physicians and clinical files. The theoretical time of the end of a wave is therefore t=+∞. The parameter τ_∗ is an under estimation of this date. It is defined as: ${\mathsf{\tau }}_{\mathrm{*}}=\mathrm{i}\mathrm{n}\mathrm{f}\left\{\mathrm{t}/\left|\mathrm{X}\mathrm{m}\left(\mathrm{t}\right)-\mathrm{N}\right|<\mathsf{ϵ}\right\}$. To compute values (K0, K1, K2) of equation (3), we first began with wave 1 using the Jacquard algorithm, (Table 1). The end of wave ${\tau }_{1*}$was estimated to occur at week 37. For wave 2, the observed frequencies were therefore corrected between the following date: beginning of the wave 2: $t_{\mathrm{0,2}}=26$ and ${\tau }_{1*}=37$ (numbers express weeks). The correction involved subtracting the frequencies of wave 1 to wave 2 of those observed during the reporting period. The corrected frequencies allowed thus to take into account the overlapped of the two subsequent waves before applying the jacquard algorithm to wave 2, (Table 1). The end of wave 2 was estimated at week 62. As the observed frequencies were underestimated compared with the mathematical model (Figure 3b), no correction was necessary for wave 3. The computed parameters for wave 3 are given in Table 1. For all the waves, the quality of fit was high with a coefficient of determination R² higher than 0.997. The computation of the total number of cases:

$X_W\left(t\right)=X_{m,1}\left(t\right)I_{[0,+\mathrm{\infty }[}+$ $X_{m,2}\left(t\right)I_{[26,+\mathrm{\infty }[}+X_{m,3}\left(t\right)I_{[53,+\mathrm{\infty }[},\forall t\ge t_0$ (4), during the three waves is shown in Figure 3b.

3.6. Parameters of Interaction between Disease and Environment

K0 parameter represents the observed speed of spread of the outbreak (Table 1). Wave 1 was significantly (p<0.05) the most virulent, followed by wave 3 and wave 2 in decreasing order. These results are in agreement with the literature. EV-71 present a reproductive number (R0) close to five, whether it was estimated to 2.5 for CV-A16 [13]. With respect to patient classification, the definition of X_M(t) (3) applied to a wave W, makes possible to build the probability for an individual to belong to a specific wave. Indeed, with a specific constant A

$$f_w\left(t\right)=\frac{1}{A}{X}_m\left(t-t_0\right),\forall t\ge t_0\mathrm{with}\mathrm{A}={\displaystyle {\int }_{t_0}^{\tau }{X}_m\left(t-t_0\right)dt}$$

(4)

fw(t) is a density probability. So

$$P\left(\raisebox{1ex}{$t_0\le u\le T$}\!\left/ \!\raisebox{-1ex}{$W$}\right.\right)={\displaystyle {\int }_{t_0}^T{f}_w\left(t-t_0\right)dt},$$

is the probability to have an inclusion time between [t0, T] when a patient belongs to a given wave W. It is therefore possible to extract patients belonging to a specific wave with a certain probability and then to express the probability for an individual to belong a specific wave. We could confirm that the severity factor and the wave factor are significantly correlated (p<0,00001). The percentage of severe cases is less important in wave 1 compared to the others. It is the result of guidelines modification. The possibility to modelize wave 2 similarly to wave 1 and wave 3 confirm the positive effect of guidelines modification. More attention was paid to severe patients. We notice that the age distribution shifted on younger patients in wave 3 compared to the others (p<0,0001). This difference could be related to the etiological agent of the disease.

4. Discussion

We question the representation of the epidemic at individual, commune, district and region scale. Figure 1b shows that the three wave can be heterogenous in terms of cases severity. Guideline’s modifications published in 2012 emphasize an important point in epidemiology: an individual is taken into account only when it is recorded by the health system. We did observe some local variation, but the guidelines were very well observed, and the health system's response was rapid, even though it took place in the middle of wave 2. The multi-scale study enabled us to show how the epidemic unfolded over the 2011-2012 period. The first wave arrived from the north (supplementary figure 1). The main focus of emergence remains the city center of Hải Phòng. Waves 2 and 3 emerged in urban areas, but further south in Lê Chân and Hai An districts. These two waves followed similar diffusion route, marked by clusters 1 and 2 of communes obtained by comparing monthly distributions of the number of cases (Figure 2). The local heterogeneity motivated the search of a general method to model over time the multi-wave epidemic observed at scale of the Hải Phòng region. The question was addressed specifically to wave 2. Beside guidelines modification, wave 2 was also characterized by a high level (close to 50%) of negative PCR, suggesting the presence of specific etiological agent.

The mathematical model developed in this work is based on the virulence/transmissibility constant, in other words, on differing dynamics of the host-virus interaction. The model is based on the principle that the key factor is not the overall number of patients but indeed the differing dynamic of host-virus interaction. Compiling the number of cases has for consequence to smooth all data and to show the epidemic as unimodal. Regardless of its true nature, the epidemic is considered associated with a single causative agent. The model developed in this work allows to breakdown the epidemic into its various components when they exist and display the epidemic as a multimodal curve. In this case each mode corresponds to a specific agent displaying a specific dynamic.

The model demonstrates that three distinct waves occurred during the 2011-2012 Hải Phòng HFMD epidemic, each one displaying a specific dynamic of expansion. According to K0 parameter, the level of host-virus interaction was different for each wave. Transmissibility was therefore a discriminating factor, each wave displaying a different speed, the highest transmissibility being associated to wave 1. Wave 2 displayed an intermediate speed whereas wave 3 was characterized by the slowest speed of expansion. This result is in agreement with R0 values which are higher for EV-71 than for CV-A6 and CV-A16 [13]. Further analysis will be performed to establish the relationship between R0 and K0.

The model developed in this work also allowed to assign the beginning and end of each wave with a highly significant fit with observed data. It was thus possible to show that the three waves partly overlapped, each wave starting before the beginning of the previous one. This indicates that the replacement of causative agent occurs during the span a given wave and not after this wave. It is therefore difficult to only explain this phenomenon by an adaptation of the human population through immune-resistance. A hypothesis could be that each causative agent is capable of infecting only a specific, susceptible part of the human population and disappears when the available naive population falls below a given threshold. However, it is in this case difficult to explain why all viral strains do not expand at the same time providing they do not target the same populations and why there is a pattern of successive waves. The relative virulence of each viral strain might therefore play a key role, in particular in asymptomatic patients where competition between strains might occur. Nevertheless, this also suggests that all strains circulate at the same time and that a given strain will predominate and expand depending on the ratio immune-resistance/transmissibility. Although we cannot clearly explain the phenomenon, the model developed in this work can thus well describe the replacement of virus along the epidemic and accurately determine the start and end of each individual component of the epidemic.

Beyond that, the model is also capable of addressing the very important issue of calculating the probability of a patient to belong to a given wave and to classify the patients into specific clusters with a high level of confidence. These clusters can be associated with a specific typology of symptoms and dynamic traits. Doing so, the model allows to estimate during the span of the epidemic if a new wave occurred, before the peak of this new wave and to identify the patients affected by this new wave comparatively with patients affected by the previous wave. This in turn open ways for clinician to determine if a new treatment strategy or a new crisis management is needed and if so, provides also the means of identifying the patients to be considered in priority. The most important aspect being that this can be done during the span of the epidemic. The model provides thus means to facilitate real time actions.

Beyond the specific case of HFMD, the model described in this work can be applied to other diseases. If the disease considered is a single phase, unimodal disease caused by a single agent, the model will bring nothing more than existing single-wave models. Similarly, if the disease is a multiphase, multimodal, disease caused by different agents but displaying the same virulence/pathogenicity traits, it can thus be considered similar to a single wave disease and here also the model will bring nothing more than existing single-wave models. However, if the disease is a multi-wave disease involving causative agents characterized by differing virulence / transmissibility traits, this model will bring very useful applications for managing this epidemic and for identifying an emerging wave associated to a potentially new strain. This work and the model presented can therefore bring very valuable support to public health in the management of multi-wave infectious diseases. Indeed, it is to our knowledge the first time such a multi-wave model with a very high fit with observed data, capable of typing patients based on clinical description and determining the emergence of a new wave has been developed. Since it is not relaying on molecular or serological test, but only on clinical parameters, the model is fast and easy to implement with no delay in response. The implementation of such model on HFMD and other multi-wave diseases will therefore bring valuable support in managing these important sanitary burdens.