A realistic two-strain model for MERS-CoV infection uncovers the high risk for epidemic propagation

Middle East Respiratory Syndrome Coronavirus (MERS-CoV) causes severe acute respiratory illness with a case fatality rate (CFR) of 35,5%. The highest number of MERS-CoV cases are from Saudi-Arabia, the major worldwide hotspot for this disease. In the absence of neither effective treatment nor a ready-to-use vaccine and with yet an incomplete understanding of its epidemiological cycle, prevention and containment measures can be derived from mathematical models of disease epidemiology. We constructed 2-strain models to predict past outbreaks in the interval 2012–2016 and derive key epidemiological information for Macca, Madina and Riyadh. We approached variability in infection through three different disease incidence functions capturing social behavior in response to an epidemic (e.g. Bilinear, BL; Non-monotone, NM; and Saturated, SAT models). The best model combination successfully anticipated the total number of MERS-CoV clinical cases for the 2015–2016 season and accurately predicted both the number of cases at the peak of seasonal incidence and the overall shape of the epidemic cycle. The evolution in the basic reproduction number (R0) warns that MERS-CoV may easily take an epidemic form. The best model correctly captures this feature, indicating a high epidemic risk (1≤R0≤2,5) in Riyadh and Macca and confirming the alleged co-circulation of more than one strain. Accurate predictions of the future MERS-CoV peak week, as well as the number of cases at the peak are now possible. These results indicate public health agencies should be aware that measures for strict containment are urgently needed before new epidemics take off in the region.


Introduction
The first case of Middle East Respiratory Syndrome Coronavirus (MERS-CoV) infection was identified in Saudi Arabia in 2012 [1][2][3][4][5]. Since then the country suffers from repeated outbreaks of MERS-CoV in different provinces (Fig 1A and 1B) [2]. There are two main routes of MERS-CoV transmission: animal-to-human and human-to-human [4][5][6]. It is suspected that dromedary camels are the source of human infections [1] but the transmission route of MERS-CoV to humans is yet not well understood ( Fig 1C). However, transmission from camels to humans is confirmed from the isolation of near-identical strains of MERS-CoV from epidemiologically coexisting camels and humans [7]. The patients might be exposed to MERS-CoV by consumption of raw camel products [8] (e.g. milk and dairy products, raw meat, etc.). Meanwhile, human-to-human transmission has been reported in society and hospital settings  with the virus being transmitted between humans during close human-tohuman contact through droplets of respiratory secretions. Potential propagation to nearby and more distant regions is also a high-risk possibility as an outbreak of MERS-CoV is likely to emerge in areas such as nearby countries in the Middle East and eastern Africa where the camel trade connects the different regions (Fig 1). Unfortunately, as for now, MERS-CoV vaccines are only at the preclinical phase [10], increasing our understanding of its epidemic potential and knowledge on the drivers of MERS-CoV variability might help to achieve better preparedness ahead of forthcoming epidemics.
The ability to predict disease outbreaks will provide a mechanism for policymakers and health-care services to respond to epidemics in a timely manner, reducing the impact and maximizing the limited resources available to be deployed [11]. The timing and severity of infectious disease outbreaks, two matters of considerable public-health relevance, are the main challenges when attempting to predict disease outbreaks [11][12][13][14]. Attempts to set up prediction frameworks for anticipating epidemics for other diseases such as dengue and influenza were pursued in the recent past with different degree of success, but clear added value [11][12][13][14][15][16]. These systems proved effective to better anticipate future outbreaks and increase our understanding on mechanisms driving disease variability. To this end, a threshold quantity that is most important to anticipate the risk of future outbreaks is the basic reproduction number (R 0 ). There are quite a few studies that estimated R 0 for the current MERS-CoV outbreak. Majumdar et al. [28] for instance, estimated R 0 for Jeddah and Riyadh using a function called Incidence Decay with Exponential Adjustments (IDEA). They found that the estimate of R 0 in Jeddah and Riyadh are in the ranges (3.5-6.7) and (2.0-2.8), respectively. A stochastic epidemiological model of MERS-CoV with zoonotic and human-to-human transmission was considered by Poletto et al. [29] to quantify the rates of generation of cases from those two transmission routes. They found that spring 2014 cases led to the increase in transmission rates, which brings R 0 to values above unity. Heishet al. [32] used a simple mathematical model to trace the temporal course of South Korea MERS-CoV outbreak. They estimated R 0 to be in the range of 7.0 to 19.3. Instead, Nishiura et al. [30] estimated the expected number of secondary cases following the importation of an index case (countries other than KSA, Qatar and UAE), using a branching process type modelling approach. They also suggested that even if R 0 is below unity, a large cluster of cases with multiple generations might occur. The aforementioned studies indicated that the basic reproduction number is above unity, although there are few studies that suggest R 0 is less than one [4,17,18]. However, at the same time, intriguing differences in observed and predicted scenarios clearly suggested other important factors might be missing in previous models built for these MERS-CoV case studies. For instance, in a recent survey for Mers-CoV in over 1300 Saudi Arabian camels, Sabiret al. [1] found that camels share three coronaviruses (CoV) species with humans. Among them, one has been dominant in the region since December 2014 and led to the human MERS-CoV outbreaks occurring in 2015. The wide species range of CoVs and their propensity to cross species boundaries suggest that more MERS-CoV uncontained outbreaks may likely occur in the future. Zumlaet al. [19] also suggested that MERS-CoV species can mutate to have increased inter-human transmissibility. Cottenet al. [37] used MERS data from the period May to September 2013. They found four Saudi Arabia MERS-CoV phylogenetic clades, with 3 clades apparently no longer contributing to current cases, therefore not appearing in the saliva of camels. They also show that the ancestors of most of the viral clades originated in Riyadh (See Cotten et al. [37]). Recently, strain variability in MERS-CoV infection was confirmed in 51 respiratory samples from 32 patients, confirming that more than one strain of human MERS-CoVis present [20]. The fact that more than one MERS-CoV strain is currently circulating in Saudi Arabia should be properly accounted for as this presence may have substantially contributed to amplify the transmission intensity producing repeated MERS-CoV changes in incidence through their periodic reintroduction into human populations [19]. Furthermore, Chan et al. [34] characterized MERS-CoV viruses from dromedaries in Saudi Arabia and Egypt and compared them with a human MERS-CoV reference strain. They isolated three dromedary strains, two from Saudi Arabia and one from Egypt. The human and dromedary MERS-CoV strains had similar viral replication competence in Vero-E6 cells and respiratory tropism in ex-vivo cultures of the human respiratory tract. They also suggested that dromedary viruses from Saudi Arabia and Egypt are probably infectious to human beings [34].

a) Disease incidence functions
We constructed three different 2-strain models for MERS-CoV that consider community human-human, hospital human-human, and passive zoonotic transmission (Fig 1C and see also S1-S18 Tables). Models derived also incorporate the effect of strain variability with strain-1 as the dominant strain ( Fig 1C and S1-S18 Tables). Variability among models is based upon three different disease incidence functions (e.g. BL, NM and SAT models [21][22][23] and see formulation in Table 1A and Fig 2).
The top panel in Fig 2 represents the BL force of infection. Here, as the number of infected individuals increases the disease transmission also increases linearly. The bottom-left panel represents the NM incidence function. Biologically this incidence function can account for "psychological effects" [38][39]. In the presence of psychological effects, initially the disease transmission rate increases rapidly when the number of infected individuals is small. However, this rate falls also rapidly in the presence of a large number of infected persons in the community. As in the case of MERS-CoV infections the CFR is about 40%, psychological effects for this infection represents the effect in the community of fear of becoming infected. This fear effect reduces the transmission rate rapidly in the presence of a high number of infections in the community. The bottom-right panel of Fig 2 represents the SAT incidence function. Biologically this incidence function represents the crowding effect in disease transmission. This crowding effect explains how the number of new infections becomes constant as the number of infected individual increases. This effect is known to reproduce the awareness effect of the disease during the course of an epidemic.

b) Super-spreaders and 1-strain vs 2-strain models
Furthermore, we investigated the potential effects of super-spreading events by incorporating an additional compartment for super-spreaders to our 2-strain model. We also consider the effect of variable zoonotic transmission by incorporating dynamics of the camel population into our 2-strain model. Information about the calculation of the epidemiological parameter values for the newly proposed 2-strain models is provided in Table 1B  respectively ( Fig 1A). Main parameters of the different models are provided in Table 1B. Parameters estimates are obtained by fitting the models to new MERS-CoV hospitalized weekly data [40]. We also estimated from data itself some unknown initial conditions for the model. Delayed Rejection Adaptive Metropolis algorithm is used here to sample the 95% confidence.The goodness of fitsto compare among 2-strain and single strain models were tested by determining the respective AIC and BIC values (See S19 Table). To further compare among predictive capabilities of these models with regard to those of their single-strain versions, predictions of raw clinical cases were attempted. Albeit some discrepancies exist, fairly accurate results were obtained despite the low numbers and highly stochastic nature of the three datasets and the nearly operational capacity of these models (Fig 3). In fact, near-operational systems usually require a much longer training period than the scarce three years employed here. Differences arising from the comparison of simulated and observed cases may also be due to the fact that we are fitting the 2-strain and single strain models to cumulative MERS-CoV cases, rather than to new cases. . The procedure for model adjustment is as follows: for approaching near-operational predictions in each province, we followed former initiatives (e.g. on influenza [14] and dengue [16]) and partitioned the whole clinical cases datasets into two parts. First part was used for calibrating the models and the remaining part (52 weeks) were left aside for out-of-fit prediction. For Riyadh, Macca and Madina prediction periods were, respectively, from weeks 105-156 (7 th June, 2015 -11 th June, 2016), weeks 65-116 (11 th July, 2015 -2 nd July, 2016), and weeks 63-114 (11 th July, 2015 -2 nd July, 2016) (Fig 3). We generated predictions for three major characteristics of the epidemiological cycle similar to previous attempts made for cholera, dengue and influenza [15,[41][42][43]

c) Estimation of R 0 , R C , and R H and prediction framework
We estimated the basic reproduction number (R 0 ), the community reproduction number (R C ), and the hospital reproduction number (R H ) for the whole period of data for Riyadh, Macca and Madina (S1-S18 Tables; 44- Afterwards, we keep on adding 4 data points to the previous data and re-estimate the parameters. These estimated parameters were then plugged into obtain values of R 0 , R C , and R H in 4 th week to 48 th week of the prediction period. Thus in intervals of 4 weeks, we obtain an estimate of R 0 , R C , and R H . Temporal evolution of R 0 , R C , and R H for the three provinces is depicted in S1-S4 Figs, respectively. Using the model that provides the best prediction for the season total cases in the three provinces (i.e., a 2-strain model with BL incidence for Macca and the 2-strain model with SAT incidence for both Madina and Riyadh), we predicted the total number of cases in the following year. Using the parameter estimates in the last prediction point (i.e. 48 th week) of the current season, we draw 1000 samples. Based on these parameter samples, we determined 1000 estimates of the season total cases in the next year [43]. We then defined a large outbreak in the forthcoming season in a particular region, as those events when cases exceeded the total number of cases in the previous season by more than one standard deviation (see Fig 5). Thus, the probability of a large outbreak (P l ) in those provinces isdetermined as the ratio of samples that exceed the total number of cases in a season divided by the total number of samples. We also considered the case when two strains of MERS-CoV may be co-circulating in the human population in Saudi Arabia. We also simulated the situation when one strain is assumed to be more active with a higher transmission rate, whilst the other is much less transmissible among individuals in the different provinces of Saudi-Arabia [37].Although we have considered two strains circulating in the community setting, we do not distinguish among strains in hospital premises. This is due to the lack of strain specific data in hospital settings. However, we are able to distinguish the contribution of infection from Community and Hospital setting by estimating the Community reproduction number (R C ) and the Hospital reproduction number (R H ). However, to account for this effect, we propose three different forces of infection (BL, NM and SAT) to model the MERS transmission in the three provinces. Saturated (SAT) functions describe "Crowding effects" in disease transmission, whereas Non-monotone (NM) functions are used to incorporate"psychological effects". As the CFR in MERS-CoV is about 40%, it is easily expected that awareness in front of a MERS-CoV situation may play an important role during an epidemic [44][45][46]. Thus a model with such a crowding effect might seem a priori more realistic in comparison to the bilinear (BL) transmission (see also Fig 2). The crowding effect in disease transmission [38,39] can also be interpreted as the behavioral change of susceptible individuals when the number of infected individuals increases (Fig 2). This fact causes a lower number of new infections when a large number of infected individuals are already present in the community. This behavioral change in the susceptible population may occur due to, for instance, public health awareness campaigns and alerts raised by local or international authorities. The psychological effect is somewhat similar to the crowding effect, but with the effect of fear being added to the awareness. The psychological effect is stronger as the transmission rate decreases more rapidly with an increasing number of infected individuals (Fig 2) [40][41][42]. Recently, Kucharski et al. [35] suggested that MERS-CoV transmission is over dispersed and hence outbreaks can include super-spreading events. In a similar study, Nishiura et al. [33] concluded that super-spreaders who visited multiple healthcare facilities drove up the epidemic by generating larger number of secondary cases. Therefore, we also consider the role of potential super-spreaders in the context of a 2-strain model with BL incidence function.

Results and discussion
Estimates of the transmission rates for the 2-strain model with BL incidence suggest that MERS-CoV transmission in the three locations is dominated by community and hospital transmission (S1 Table to S9 Table). The former statement is in good agreement with a previous study that suggested that MERS-CoV infections are essentially produced through both hospital and community based human-human transmission [2]. This fact is well reflected in the estimates of transmission rates obtained from the other two 2-strain models (e.g. NM and SAT incidence functions) for the three locations (S1-S9 Tables). Estimates of the parameter θ (a measure of variability within the two strains) in Riyadh and Macca (see S1-S3 Tables) suggest that both strains are currently active in these provinces (albeit strain 2 contributes with only 3% -50% to community human-human transmission). In Madina instead, a majority of the models (BL and NM) point to only one active strain,with strain 2 contributing less than 1% to community human-human transmission, see S1-S9 Tables.
For each particular region in this study, two strain models are compared among themselves and also with regard to their single-strain counterparts in order to determine the best model (i.e. single, two-strain or a combination of other two-strain models). Results show that in the three locations, 2-strain models provide a better fitting to cumulative cases in comparison to their single strain counterparts (Fig 2). AIC and BIC values for the 2-strain models suggest that the best model is region specific (best model for Riyadh is a 2-strain model with BL incidence, whereas in Macca and Madina the best models are both the model with a NM incidence function and the one incorporating super-spreaders, respectively; see S19 Table). We also compare our 2-strain models with a 2-strain model with variable Zoonotic transmission (SI Eq. A2). Comparing AIC and BIC values, we found that the models with variable Zoonotic transmission (SI Eq. A2) did not improve all the previous results (see S19 Table). Therefore, the best model among all 2-strain models cannot be determined solely from their goodness-offit comparison. For this reason, we additionally compared all the two-strain models on the basis of their respective prediction skills. Forecasts for the three models in Riyadh suggest that, a majority of the times, the 2-strain model with SAT incidence can better predict the three targets in comparison to the other two competing 2-strain models (e.g. note that for SAT the average of the mean absolute error, MAE, for peak week is 11.6, whereas for peak maximum is 24.83 and for season totals 55.62) (Table 2A). For comparison, we also provide predictions in the province of Riyadh of the three targets using the model with super-spreaders. However, prediction with super-spreaders did not improve at all the former results obtained for Riyadh (see Table 2A). For Madina, both the peak week and season totals are better predicted by the same SAT 2-strain model (i.e. MAE average for peak week being 14.5 and 9.33 for season totals; see also S22 Table). However, when predicting the peak maximum in Madina, clearly the best model is the 2-strain model with the NM incidence function (i.e. MAE average for peak maximum being 2.62; see also S22 Table). For Macca, except for the peak week, the other two targets, namely peak maximum and season totals are better predicted with the 2-strain model with BL incidence (i.e. MAE average of 7.08 for peak maximum and of 63.42 for season totals, respectively; see S22 Table). For peak week in Macca, again a better prediction is achieved by the 2-strain model with SAT incidence (i.e. MAE average for peak week being 23.9, see S22 Table).
When comparing the three two strain models to their single-strain counterparts, composites of the mean absolute error (MAE) in Riyadh (Table 2A, and S23 Table) suggest that 2-strain models always outperform the single strain versions in predicting the three targets (i.e. peak week, peak maximum, and season totals). Averages of MAE for peak week, peak maximum and for season totals are 11.6, 24.83 and 55.62, respectively. This fact reinforces our earlier inference, confirming that there is more than one strain currently active in Riyadh. Among the single-strain models in Riyadh, the model with SAT incidence is found to have a better predictive capacity for the three targets in comparison to the other two single-strain models (S23 Table). In Macca instead, the single-strain model with SAT incidence provides better predictions of the peak week, the peak maximum and seasontotals in comparison to all the other models (see for instance MAE values in S22 Table; and values obtained for 1-strain and 2-strain models in S22 Table). Conversely, in Madina, both 1-strain and 2-strain models with SAT incidence provide a similar performance in predicting the three targets (e.g. MAE for peak week 16.7 in comparison to 14.5, for peak maximum of 3.33 compared to 3.6 and for season totals of 5.87 against 9.33; see S22 and S23 Tables). The latter may likely represent the best choice of models to form the basis of a future early-warning system for MERS-CoV prediction in the region.
The best model, selected on the basis of its capacity to derive skillful predictions, is the model with the SAT incidence function. We used this model to estimate the basic reproduction number, R 0 , the community reproduction number, R C , and the hospital reproduction number, R H , for the whole period of data available in those three provinces (Table 2B and S20  Table). In all three provinces, R 0 is estimated to be always greater than unity. In all three provinces, R H is found out to be larger than R C . This implies, MERS-CoV transmission triggered from the hospital setting (see Table 2B). These estimates agree with some previous values reported in the literature [4,[17][18], while being a result well supported by other R 0 estimates [29][30][31]33]. However, previous modeling attempts to model Saudi Arabia MERS-CoV clinical incidence used only a 1-strain model with BL transmission [4]. In our case, versions of R 0 were also estimated from our BL 1-strain model in those same three provinces (see S21 Table) and they are in good agreement with the previous values provided by Chowell et al.[4]. At the same time, though, 1-strain model show less predictive capacity than their similar two-strain counterparts. To further verify the robustness of these estimates, the temporal evolution of R 0 , R C , and R H are displayed for different time intervals in the three provinces (S1-S4 Figs). Considering the best model configurations, we additionally computed the temporal evolution of R 0 , R C , and R H in the three provinces. Temporal changes in R 0 (S1 Fig) indicate that in most of the predicted weeks, R 0 stays well above the epidemic threshold (R 0 = 1). This fact is well established from the temporal evolution of R C , and R H in those three provinces (S1 and S2 Figs).

Table 2. (A)Average predictions [Simple average of Mean Absolute Errors (MAE)] over all forecasts of the three 2-strain models and their different combinations for the Riyadh province. (B)Estimated values of the Basic reproduction number (R 0 ), the Hospital reproduction number (R H )
, and the Community reproduction number (R C )for the three provinces of Saudi Arabia using best model (two strain). The best two strain MERS model is with saturated incidence. The data are given as Mean (95% CI). � The MAE for the three models and their combinations during point prediction of peak week and peak incidence were calculated up-to the peak of the prediction season (11th Week) of the Riyadh province. Seasonal forecasting of MERS-CoV epidemics As the basis for an operational EWS for the region, we predicted the total number of cases in the out-of-fit interval covering July 2016 to July 2017. We used the best individual model for the three provinces among all 2-strain and single strain models developed. Figs 6 and 7 and Table 3 clearly indicate the 2016-2017 season might be ripe for a larger outbreak in Macca and Riyadh. Instead, in Madina the likelihood of suffering a larger outbreak is very low (Fig 6 and Table 3). Fig 7 displays the single-strain model fitting to new MERS-CoV cases. Here M1 refers to the single strain model with BL force of infection, M2 to the single strain model with nonmonotone incidence, and M3 to the single strain model with saturated incidence. Solid black curve represents model solution and yellow region denotes 95% confidence interval for predictions. Simulations for both Macca and Riyadh reflect quite appropriately the dynamics displayed in observations (Fig 7). Fig 8 shows All-season's predictions for the three consecutive years from July 6, 2013 to June 28, 2016, with each interval of 52 weeks fitted shown in panels S1-S3.
According to WHO reports, 249 MERS-CoV cases including 75 deaths (CFR 30%) were reported from Saudi Arabia between July 2016 and July 2017 [44]. In the aforementioned period, at least 108 new cases [44] were reported from Riyadh province and at least 5 cases were reported from Madina. These values indicate that a larger epidemic did not occur in the last season (2016-2017), the risk for a larger one in the coming seasons still remains high.

Conclusions
Up to July 2018, 2237 new cases were reported, with 1861 only in KSA and 793 deaths [45] (CFR 35.5%). This is an alarming situation as previous predictions on MERS-CoV had instead suggested that MERS-CoV might not sustain as an epidemic in the Arabian Peninsula. The WHO report [44] suggests that MERS-CoV is still a relatively rare disease about which the medical personnel in health-care facilities have low awareness. Globally, MERS-CoV Seasonal forecasting of MERS-CoV epidemics awareness is limited and because symptoms of MERS-CoV infection are non-specific, initial cases can be sometimes easily missed. With improved compliance in infection prevention and control, namely by stricter adherence to the standard precautions at all times, human-tohuman transmission in health-care facilities can be reduced and even possibly eliminated with additional use of transmission-based precautions. In that regard, predictive mathematical models can help strengthen our understanding of both MERS-CoV transmission and control.
In this study, we addressed the capacity of predictive mathematical models based on twostrain MERS-CoV configurations having different transmission functions. The models differed from each other in their force of infection and in how they cope with heterogeneity in transmission. Estimates of transmission rates suggest that community and hospital transmission are dominant in the case of 2-strain models in Riyadh, Macca and Madina. The majority of the 2-strain models suggest that MERS-CoV transmission is dominated by community and hospital human-human transmissions, a fact that reflects the actual transmission scenario in Saudi Arabia [2][3][4]. Estimates of the parameter that measures transmission diversity between the two strains in the three provinces suggest that two strains are only active in Riyadh. This opposite trend in Riyadh in comparison to the other two provinces may be due to the fact that Riyadh is  Seasonal forecasting of MERS-CoV epidemics the capital city with good large health care facilities and a majority of the MERS-CoV patients in Saudi Arabia come to Riyadh for treatment [2]. These patients may therefore carry different MERS strains, ultimately leading to multiple strains being presently co-circulating in Riyadh. Similarly, Cotten et al. [37] found that ancestors of most of the viral clades originated from Riyadh. We compared among the 2-strain models according to their predictive performance with regard to three targets (i.e. peak week, peak maximum and season totals). Our results suggest that among the three 2-strain models, the model with SAT incidence provides consistently skillful predictions and may be used to date as the best predictive model for MERS-CoV in Riyadh. Riyadh iswhere most of the MERS-CoV cases occur, while for the provinces of Macca and Madina, with lower reported MERS cases, it is difficult to determine the best model among the three 2-strain models. This fact justifies our earlier finding that in Riyadh two different strains are currently active and therefore the performance of the 2-strain models is better there. As per our results for Macca and Madina, only one dominant strain is active in those provinces. Therefore, predictions based on single strain models are there more appropriate. Our results also suggest that among the single strain models, those with SAT incidence always accurately predict the three targets for these two provinces. Thus, a dynamical MERS model considering this crowding effect is the most appropriate configuration to cope with the nature of MERS-CoV transmission.
We estimated R 0 using the best 2-strain model in Riyadh and estimates are in good agreement with the findings of Majumdar et al. [28].The finding that R 0 is most of the time above 1 . One year appears in each panel S1 to S3(S1 = S2 = S3 = 52 weeks). Line and circle line refers to observations and 95% confidence interval for cumulative predictions denoted by the yellow region.
https://doi.org/10.1371/journal.pntd.0008065.g008 (S1 Fig) is well supported by some previous estimates in literature [28][29][30]32]. Lower contribution of community transmission in R 0 (See Table 2B) in Riyadh and Macca suggests that MERS-CoV transmission is triggered from hospital settings in those provinces. Most interestingly, in some forecasted periods, R 0 attains large values, a fact that denotes that a rapid propagation among the susceptible population is indeed possible. More worrisome is the range of values into which R 0 moves, most of the time above 1 and below 2,5, a fact that makes it a dangerous infection in terms of silent and constant potential population spread. Community reproduction number R C well above unity (see S4 Fig) in most of the predictive weeks indicate that in a near future a large outbreak may be possible in those provinces. Out-of-fit predictions for the next season totals suggest that there is a high possibility of larger outbreaks in Macca and Riyadh. However, our results instead indicate that there is a very low possibility of larger outbreaks in Madina. The fact that this outbreak did not happen in 2017-2018 does not preclude what may occur in the forthcoming seasons. Under such a scenario, authorities and international health agencies should prepare and actively work towards the prompt implementation of cheap albeit efficient computational platforms ready to assist in the simulation of how a potential outbreak might evolve in the region. More so, given the high probability that another large MERS-CoV outbreak occurs in the Arabian Peninsula or nearby countries. Migration may also play a major role for increased transmission in the provinces of Saudi Arabia and this feature should be properly accounted for in future model configurations. In summary, our findings suggest that in a majority of provinces a single MERS-CoV strain is currently active, conversely to the situation in Riyadh. However, in the near future, it is also possible that more general MERS-CoV transmission occurs from multiple strains in other provinces of Saudi Arabia.
Supporting information S1