Modeling , analysis and simulations of MERS outbreak in Saudi Arabia

This work describes a continuous differential equations model for the dynamics of Middle Eastern Respiratory Syndrome (MERS) and provides its computer simulations. It is a continuation of our previous paper Al-Asuoad et al. (Biomath 5, 2016) and it extends the simulations results provided there, which were restricted to the city of Riyadh, to the whole of Saudi Arabia. In addition, it updates the results for the city of Riyadh itself. Using an optimization procedure, the system coefficients were obtained from published data, and the model allows for the prediction of possible scenarios for the transmission and spread of the disease in the near future. This, in turn, allows for the application of possible disease control measures. The model is found to be very sensitive to the daily effective contact parameter, and the presented simulations indicate that the system is very close to the bifurcation of the stability of the Disease Free Equilibrium (DFE) and appearance of the Endemic Equilibrium (EE), which indicates that the disease will not decay substantially in the near future. Finally, we establish the stability of the DFE using only the stability number Rc, which simplifies and improves one of the main theoretical results in the previous paper. Keywords-MERS model; stability of DFE and EE; simulations; sensitivity analysis;


I. INTRODUCTION
This work uses the mathematical model constructed in [2] to study the dynamics of the Middle East Respiratory Syndrome (MERS) in Saudi Arabia.It also expands the study that was performed there of the disease in the city of Riyadh, since new data became available since the publishing of the paper.The aim of this work is to provide the health care community and related authorities with a predictive tool that allows to assess various MERS scenarios and the effectiveness of various intervention practices.
MERS is a new respiratory disease caused by the newly discovered Middle East Respiratory Syndrome Coronavirus (MERS-CoV).The first case of the disease was reported in Saudi Arabia in June 2012, when a 60-year-old man died of progressive respiratory and renal failure 11 days after hospital admission.He had a 7-day history of fever, cough, expectoration and shortness of breath [26], [30].
In September 2012 a case of a 49-year old man from Qatar with pneumonia and kidney failure, who was treated in an intensive care unit at a London hospital, was reported.He had a history of travel to Saudi Arabia.Further laboratory tests revealed a positive MERS-CoV infection [6].Retrospectively, the infection was found in stored respiratory and serum samples on two deceased patients from Jordan, where in April 2012 an outbreak of acute respiratory illness occurred in a public hospital [17].
In 2003, a previously unknown coronavirus, the Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV), caused a global outbreak of pneumonia that resulted in approximately 800 deaths [22].MERS-CoV that has effects similar to those of SARS-CoV, is classified as a coronavirus, which is a family of single-stranded RNA viruses.This family includes viruses that cause mild illness such as common cold as well as severe illness such as SARS in humans.MERS-CoV is a beta coronavirus which has not been identified in humans before 2012 and is different from any coronaviruses (including SARS-CoV) that have been found in humans or animals [10], [31].
Within a year from its discovery, a total of 130 MERS-CoV cases were identified, 58 of which died, which means that the case fatality rate is 45%, much higher than SARS-CoV, which has a case fatality of 15% [24], [27].Up to date (August 31, 2017) 2067 confirmed cases of MERS-CoV have been reported worldwide, out of these 1679 were reported from Saudi Arabia where the case fatality rate has been 40.6%[28].The infection has been a global threat due to continuous outbreaks in the Arabian Peninsula and international spread to 26 countries including Qatar, Jordan, United Arab Emirates, United Kingdom, the Philippines, United States and other countries, by infected travelers [14].In 2015, a large outbreak happened in South Korea, which was the first outbreak outside the Arabian Peninsula [19], [23], and 186 people were infected, 38 of whom died.
After intensive search, camels were found to have a high rate of anti-MERS-CoV antibodies, which indicates that they were infected with the virus.Then, definite evidence of camel-to-human transmission of the virus has been reported recently [5], [29].Moreover, there is clear evidence that the infection is transmitted from person to person upon close contact, including from patients to healthcare workers [4], [12].
The incubation period from exposure to the development of clinical disease is from five to 14 days.MERS-CoV is typically characterized by cough, fever, sore throat, chills, myalgia and shortness of breath [11], [13].One-third of the patients had also gastrointestinal symptoms such as vomiting and diarrhea.The common complications of the MERS-CoV infection include pneumonia, acute respiratory distress syndrome and respiratory failure.Although it is known that asymptomatic infection occurs, the percentage of patients who have it is unknown, yet, [3], [7].
No specific treatment is available for the MERS-CoV infection.Currently, the management of the disease is done by supportive therapy that minimizes the symptoms.Some patients require mechanical ventilation or extra-corporal membrane oxygenation.Since no vaccine exists for MERS [9], [18], once a case is identified, the individual and those connected to them are being isolated to minimize the spread of the disease.
Part of the content of this article can be found in the recent Doctoral Dissertation [1] where additional information can also be found.
To help assessing the threat of the spread of MERS, we constructed a mathematical model as a tool to predict possible future scenarios and the effectiveness of various intervention procedures.The model is in the form of a coupled system of five nonlinear ordinary differential equations (ODEs) for the susceptible, asymptomatic, clinically symptomatic, isolated and recovered populations.It is of a rather standard MSEIR type (see, e.g., [8], [15], [16], [21] and the many references therein).The novelty in this work lies in the theoretical proof of the global stability of the Endemic Equilibrium (when it exists), and simulations for the whole of Saudi Arabia.
First, we mathematically analyze the model and provide a new proof of the global stability of disease free equilibrium (DFE) and of the endemic equilibrium (EE), when it exists (and then the DFE becomes unstable).Moreover, the new proof of the global stability of the EE is simpler and more elegant than the one in Proposition 4 in [2], since it uses only the effective reproduction number or the stability control parameter R c .
Second, the model is used for simulations of the overall outbreak in Saudi Arabia, and it also extends the study of MERS done in [2] of the city of Riyadh, as more data has been collected since its publication.Indeed, currently there exists data that spans 1550 days, [25].For the whole of Saudi Arabia, the model parameters were fitted to the data using the first 865 days, then runs for 1690 days, until Nov. 4, 2020, were performed, allowing for the prediction of the disease spread in the next three years.In the previous paper, we fitted our model to the daily reported cumulative cases of MERS data for Riyadh for the period from Nov. 4, 2013 to March 17, 2016 (865 days).Here, we fitted the model to the data from Nov. 4, 2013 to July 11, 2017 (1346 days).The model was found to be very sensitive to the scaled contact parameter that is directly related to the number of individuals a person is in contact with each day.Nevertheless, the simulations provide a very good fit with what has been observed and are similar in their predictions of the near future, say the next two years.
The rest of the paper is structured as follows.The model is described in Section II, where its compartmental structure and flow chart are also provided.The stability analysis of the DFE and EE is done in Section III, where the local and global stability of the DFE are studied.Our new results on the global stability of the EE are summarized in Proposition 4 and proved using the Lyapunov method and LaSalle's principle.The description of the simulation results can be found in Section IV.First the optimized parameters for the baseline for Saudi Arabia are presented and the simulation results depicted.Then, the extended study of Riyadh is described.The sensitivity of the model to the contact parameter β is done in Section V, which is one of the main characteristics of the model.In Section VI we depict graphically the errors, i.e., the difference between the data and model predictions.The paper concludes in Section VII where the results are summarized and some unresolved issues indicated.

II. THE MODEL
We use the basic model of MERS dynamics developed in [2], where full details of the model and its underlying assumptions can be found so, the description of the model here follows very closely the one in [2].The model represents the disease dynamics of five populations of individuals: susceptible S(t), asymptomatic E(t), symptomatic I(t), isolated J(t), and recovered R(t), where the time t is measured in days.Also, N (t) denotes the total population at time t and is given by The model flow diagram is depicted in Fig. 1.The mathematical model for the MERS disease (the basic model in [2]), which consists of five rate equations for the dynamics of S, E, I, J and R, is as follows.
Model 1: Find five functions (S, E, I, J, R), defined on [0, T ] and values in R + ∪ {0}, such that for 0 < t ≤ T , the following hold: together with the initial conditions Here, we denote by S 0 > 0 the initial population before the disease outbreak, and E 0 , I 0 , J 0 and R 0 are nonnegative populations that satisfy (1) at t = 0.It is appropriate, within the context of Saudi Arabia to assume that that initially S 0 = N (0) > 0, and the others vanish, meaning that at first there are only susceptibles in the population.However, for the sake of generality, we assume that the initial populations are nonnegative.
The rate of change of the susceptible population S(t) is given in equation (2), where P represents the recruitment rate and is assumed to be a constant.We denote by β (1/day) the effective contact rate.The dimensionless parameters E and J are the transmission coefficients from asymptomatic and symptomatic individuals, respectively.Thus, the second term on the righthand side of equation (2) describes the rate at which the susceptibles become infected with the virus as a result of contact with asymptomatic, infected, and isolated individuals.The population's natural death rate is µ (1/day).Equation (3) describes the rate of change of the asymptomatic population E(t).These individuals carry the virus but have not yet developed clinical symptoms of MERS, which means that they can infect susceptibles unintentionally.Following [2], k (1/day) is the rate of development of clinical symptoms in asymptomatic population.The last term on the right-hand side of (3) describes both the mortality rate due to development of clinical symptoms at rate k and the natural mortality rate.We turn now to equations (4) and (5).It is assumed that infectives are isolated at rate γ (1/day).The parameters σ 1 , σ 2 (1/day) denote the recovery rate of symptomatic and isolated populations, while d 1 , d 2 (1/day) are the disease-induced death for symptomatic and isolated populations, respectively.The first term on the right-hand side of the equation ( 4) represents the asymptomatic individuals who developed clinical symptoms and become infected, while the first term on the right-hand side of the equation ( 5) represents the isolatedinfected individuals.Finally, the rate of change of the recovered population R(t) is given in equation (6) where the first and second terms on the righthand side represent the recover-infected and the recover-isolated individuals, respectively.We note that since there is no data, yet, about possible reinfection of the Recovered, we assume that they are permanently immune.
We recall that µ denotes the natural death rate.Thus, if a person has a a life expectancy of 80 years, then the natural death rate µ is 0.000034 per day.It was assumed that in the absence of disease, the total Saudi population was N = P µ = 32 million for P = 1088 people and µ = 0.000034 per day and the total population of Riyadh was N = P µ = 5 million for P = 170 people and µ = 0.000034 per day.
A full description of the variables, parameters, and the parameters' values considered in the model can be found in Table (I).
Finally, the cumulative cases of MERS up to time t, were obtained from the expression with the initial value CT (0) = 0, while the cumulative recovered from the disease up to time t, were obtained from Biomath 7 (2017), 1802277, http://dx.doi.org/10.11145/j.biomath.2018.02.277 with the initial value CR(0) = 0. Similarly, the cumulative deaths induced by the disease up to time t, were obtained from with the initial value D(0) = 0.
We note, for the sake of completeness, that the following results, in addition to those mentioned and improved in Section III, were established in [2]: the existence and uniqueness, positivity and boundedness of the solutions.Indeed, it was found that all the trajectories of the system lie in the set which is invariant and compact.

III. STABILITY ANALYSIS
We first analyze the stability of the Disease-free Equilibrium (DFE) and of the Endemic Equilibrium (EE), both locally and then globally.However, we note that the local stability of the EE has already been done in [2].These results improve considerably the results there, and also simplify them as they show that the effective reproduction or control number R c controls the stability, and this closes a gap described there.Thus, there is no need for the basic stability number R 0 that was introduced there.We have, Here, λ 1 is the largest eigenvalue of the Jacobian matrix J(P 0 ) given shortly, and

A. Local stability of the DFE
We begin with the local stability.It is straightforward to see that the DFE is given by P 0 = (S 0 , 0, 0, 0, 0), and S 0 = P/µ.
In our previous paper ( [2]), the local stability of the P 0 was proved in term of R 0 .Here, we used the Routh-Hurwitz criterion (see e.g., [21]) to prove the following stability result using R c .Our local result is the following.
Proposition 2: The disease-free equilibrium of the model is locally asymptotically stable when R c < 1 and is unstable when R c > 1.
The characteristic equation is where Now, since (λ + µ) 2 = 0, there are two equal and negative eigenvalues, λ 1,2 = −µ.The remaining three eigenvalues are determined from the cubic equation It follows from the Routh-Hurwitz criterion ( [21]) that the solutions of this equation have negative real parts when A > 0, B > 0, C > 0, and AB > Biomath 7 (2017), 1802277, http://dx.doi.org/10.11145/j.biomath.2018.02.277 C. Clearly, A, B, C > 0 when R c < 1.It remains to show that AB > C. We have, ).Thus, AB > C, when R c < 1, and it follow from the Routh-Hurwitz criterion that all the eigenvalues have negative real parts in this case.We conclude that P 0 is locally asymptotically stable.

B. Global stability of the disease-free equilibrium
The global stability of the disease-free equilibrium is shown in the following proposition, based on the construction of an appropriate Lyapunov function and the use of LaSalle's invariance principle.
Proposition 3: The disease-free equilibrium of the model is globally asymptotically stable in R 5 + when R c ≤ 1.
Proof: To show the global stability of the disease-free equilibrium P 0 , we construct the following Lyapunov function in which we only considered the variables representing the infected components of the model, where, Next, we let where Calculating the derivative of L along the solution (E(t), I(t), J(t)) of the system (3)-( 5), we obtain Therefore, since all the parameters are nonnegative, dL dt ≤ 0 when R c ≤ 1.We note that dL/dt = 0 if and only if E = I = J = 0 i.e., it vanishes only at the disease-free equilibrium.Therefore, if we let then the largest compact and invariant set in Γ is the singleton {P 0 }.By LaSalle's invariance principle ( [20]), every solution of the equations ( 2)-( 6), with initial conditions in Ω N , approaches P 0 as t → ∞, whenever R c ≤ 1.This completes the proof.

C. Global stability of the endemic equilibrium
The global asymptotic stability of the endemic equilibrium when it exists, is also proved by constructing an appropriate Lyapunov function and using LaSalle's invariance principle.We note that it follows from the local stability analysis (see [2][Proposition 5]) that the EE exists and is unique only when R c > 1.Therefore, we deal with the case R c > 1.
Theorem 4: The the endemic equilibrium P * = (S * , E * , I * , J * , R * ) exists and is globally asymptotically stable when R c > 1, and does not exist when R c < 1.
Proof: To show the global stability of the endemic equilibrium P * , we consider change of variables and construct the following Lyapunov function L, where, We note that L(0, 0, 0, 0, 0) = 0 and L(W 1 , W 2 , W 3 , W 4 , W 5 ) is positive.Calculating the derivative of L about the system (2)-( 6), we obtain Then, using the equations, we obtain It follows that Hence, each term of dL/dt < 0, and thus the largest invariant set at which dL/dt = 0 is the equilibrium point P * and it follows from LaSalle's invariant principle [20] that P * is globally asymptotically stable.

IV. NUMERICAL SIMULATIONS
We turn to describe the numerical simulations of the model.We used the same numerical algorithm that was developed and implemented in MAPLE in [2].Then, we run extensive numerical simulations, using the values of the parameters given in Table I for the baseline simulations.A number of other sets of parameters were also used, as explained below.The simulations were run for both the city of Riyadh and the whole of Saudi Arabia.Those for Riyadh were an extension of the simulation in [2], Biomath 7 (2017), 1802277, http://dx.doi.org/10.11145/j.biomath.2018.02.277 since new cases were found since the last day of simulation reported there, and the additional data has been incorporated into the simulations below.The simulations of Saudi Arabia were new and motivated by the fact that the model predictions in [2] were very close to what subsequently has been observed in the field.
The data currently available for the outbreaks of MERS in Riyadh and in Saudi Arabia was from November 4, 2013 till February 1, 2018, a total of 1550 days (∼ 51 months).The parameters P and µ, which were not associated with the disease, were readily available for both places.The other model parameters were obtained by fitting the numerical solutions to the data of the first 865 days of the disease breakout using the optimization routine lsqcurvefit in MATLAB.This generated the baseline case parameters that are provided in Table I.
Then, we run the simulations for a longer time and observed the model predictions in the following 540 days for which data were available but not used in the parameters fitting.This provided an insight into how well the model predicted the disease dynamics.We would like to point out here, as was noted in [2], that the additional data was found to fit very well into the model, without any need to change the previously fitted parameters, and we describe these results in detail below.
We first present the baseline simulations for the whole country, and this is completely a novel addition to the literature.Then, we study the disease spread in the city of Riyadh, where additional information was provided.Finally, we perform a reduced sensitivity analysis for the model with respect to the scaled effective contact number β, which shows that the simulation results are extremely sensitive to its value.We discuss it in Section V.
A note on the optimization for the model parameters.The optimization program found a number of local minima that provided, for the cases of Saudi Arabia and Riyadh, results in which R c had a value close to one, both below and above one.We chose the baseline simulations in both cases to be those with R c < 1, since these lead to a slightly better fit.However, below we depict simulation results for Saudi Arabia with either R c = 0.99704 for the baseline case, in which the DFE is stable and attracting, or R c = 1.004, in which the EE is asymptotically stable and the DFE unstable.Similarly, for the city of Riyadh, we used the baseline case with R c = 0.9928, which is related to a stable and attracting DFE, or R c = 1.0045 that has an unstable DFE and stable and attracting EE.These are directly related to the sensitivity of the model to β and as we explain below, it was found that a change of 0.3% leads to the change in the stability, hence in the disease dynamics.We note in the cases when R c < 1, when there in EE, the decay to the disease free equilibrium to the DFE is slow, over a period of more than a hundred years, and since our interest in this work is only the next few years, we do not depict the long time results.

A. Baseline simulations -Saudi Arabia
In this subsection, we describe the baseline simulations for the whole of Saudi Arabia with R c = 0.99704.Then, for the sake of completeness, below we describe another set of simulations with R c > 1.Both agree well with the data in the short term (1550 days ≈51 months), but differ in long term behavior, as was to be expected, since the baseline is associated with R c < 1 in which there is no EE, while the second set with R c > 1 is associated with stable and attracting EE.However, both values are very close to 1.
The daily reported new cases of MERS were obtained from the Saudi Arabian Ministry of Health website [25].More specifically, we considered the period of 1550 days from Nov. 4, 2013 to February 1, 2018.A nonlinear least square fit, using lsqcurvefit -a Matlab function contained in the optimization toolbox-was performed to obtain the model parameter values.As was noted above, we fitted the basic model parameters of ( 2)- (6) to the data from Nov. 4, 2013 until Mar.17, 2016 (a period 865 days) using reasonable initial guesses for the parameter values and obtained better estimates of the same parameters from the  optimization fit, which are given in Table I.The results of model fitting, which is the baseline, are depicted in Fig. 2.  I.
We next describe the simulations of the MERS model, equations ( 2)-( 6) with the initial condi-tions S(0) = 31, 999, 990, E(0) = 0, I(0) = 10, J(0) = 0, R(0) = 0.This choice of these initial conditions was made based on the data or the lack of it on Nov. 4, 2013 when it became available and when the simulations start.The results of the numerical simulation, depicted in Fig. 3, show a very good agreement between the model predictions-smooth colored curves-and the observed data -red dots (red curves on this scale).We emphasize again that the curve fitting was done on the first 865 days and the next 683 days are the model predictions, and they agree with the field date very well, indeed.Then, in the figure we depict the model prediction for another three years, or 1006 days, until Nov. 4, 2020, which means cumulative results for 2555 days of simulations.
In the simulations, Fig. 3, the cumulative number of: infected cases is depicted in the top (T), the recovered in the middle (M), and the deaths on the bottom (B).The model predicts, that if the epidemic continues its current trajectory, by Nov. 4, 2020 (another 33 months), there will be about 2200 new cases (M), the cumulative recovered will be about 1449 (M), and the cumulative deaths will be around 760, (B).
We note that there is an under-reporting issue with the cumulative number of recovered (M) for  the first 183 days, since the data was not available, so the number was set as zero and this explains why the whole red graph is below the blue curve.However, by raising the red dots curve to agree with the blue curve on day 183 led to a very good fit on the cumulative recovered, too.
Although in this case the DFE is stable and attracting, and the disease will die out in about 20 years, we did not show the long term behavior since at this stage it seems not to be very relevant.
It is seen that if MERS continues in the current trajectory, in the next three years one can expect another 548 cases or so in the whole country.Although the number is not large relatively to the size of the population of the whole country, the possible epidemic-like spread of the disease must be taken into account by the authorities.We return to this point in the conclusions section.
Next, as was noted above, running the optimization program with different initial conditions yielded a number of sets of values, related to local minima of the optimization function.So for the sake of completeness, we present simulation results with somewhat different parameters, provided in Table II in which case R c = 1.004.Thus, the EE is stable and attracting, and the disease cannot be eradicated.
We note that in the case when the DFE is asymptotically stable we had β = 0.1818 and here we have β = 0.1284, which is very interesting and this point is discussed further in Section 5, when we study the sensitivity with respect to β.Since in this case R c = 1.004, the endemic equilibrium exists and is asymptotically stable.Indeed, we found that the EE values P * = (S * , E * , I * , J * , R * ) were (31, 772, 111; 27; 51; 13; 113, 994).
Note that the population of the country was taken as |P * | = 31, 886, 195.
What these numbers mean is that when the disease is close to the EE, one would have on each day about 27 asymptomatics, 51 people with the disease symptoms, 13 isolated, and 113, 994 had just recovered.Clearly these numbers, if MERS would takes such a turn, pose significant challenges to the authorities and the whole society.
The eigenvalues corresponding to the Jacobian matrix J(P * ) were found to be − 0.000034, −0.399388, −0.225192 − 0.000017 ± 0.0001127i, indicating that the EE is locally stable and attracting, as was claimed above.
We solved the system with the same initial conditions as above.The simulations results are depicted in Fig. 4, where the run was for 2555 days (∼ 84 months).We note that the endemic equilibrium is approach in about 80 years.
The cumulative infected cases of MERS are depicted in the upper left (T), the cumulative number recovered in the upper right (M), and the cumulative number of deaths on the bottom (B).If the epidemic switches to such a trajectory, by the year 2020 (another 33 months), the model predicts a bit over 4000 new cases, the cumulative death will be around 2000, and there will be about 2000 recovered.
It seems, comparing the two simulations, that the fit for R c = 0.99704 is better than the one with R c = 1.004, however, visually it is also related to the vertical scales in the figures.The difference in the fit is actually quite small, although the difference in the predictions is larges.
The simulations in Saudi Arabia show that eventually the disease will either disappear, or stabilize.At this stage it is impossible do decide which is the case based on the model predictions.However, in either case, in the next few years MERS will be spreading, possibly between 620 and 1240 deaths Biomath 7 (2017), 1802277, http://dx.doi.org/10.11145/j.biomath.2018.02.277 per year.This we would like the authorities to be aware of it.

B. Baseline simulations -Riyadh
We turn to the simulations of Riyadh, and as was noted above, these include more data that became available since our paper [2].Moreover, the parameters are somewhat different than those in the article.Indeed, as was the case with Saudi Arabia, we present simulations with two sets of parameters obtained by the optimization subroutine.The first, which seems to be a better fit with the field data has R c = 0.9928, while the other one, whose fit is slightly worse has R c = 1.0045.Therefore, the first simulations are in the case when MERS will eventually disappear, while in the second case the endemic equilibrium exists and the disease will persist.
In the first case, the values of the parameters obtained from the fit for R c = 0.9928 are given in Table I.We solved the MERS model ( 2)-( 6) with the initial conditions S(0) = 4, 999, 990, E(0) = 0, I(0) = 10, J(0) = 0, R(0) = 0.The results of numerical simulation, depicted in Fig. 5, seem to agree well with the observed data.Again, we note that the fit was found using the field data reported in [25] for the first 865 days, and the very good agreement in the following 683 days just supports the model predictions.
The cumulative infected cases of MERS are depicted on the upper left (T) of Fig. 5, the cumulative number of recovered in the upper right (M), and the cumulative number of deaths on the bottom (B).If the epidemic continues at its current trajectory, by 2020 (another 33 months), the model predicts about 840 new cases, the cumulative death will be around 260, and there will be about 580 recovered.
As was done above, we now present another simulation in the case when R c = 1.0045.The values of the parameters obtained from the optimization fit are given in Table III.
We solved the MERS model with the same initial conditions and simulation results are depicted in Fig. 6.The model predicts that if the epidemic would proceed on this trajectory, by the year 2020 (another 33 months), there would be about 1440  new cases, the cumulative death will be around 740, and there will be about 680 recovered.
When comparing with the predictions of the case with R c < 1 above, it is found that there would be 480 more deaths in the next three years.
We conclude that wether the DFE is stable and attracting and eventually the disease will disappear or the EE is stable and attracting and the disease will be active for along time, in the near future, say the next three years, one can expect at least a hundred to two hundred and fifty deaths per year in the city of Riyadh.However, these scenarios depend crucially on the effective contact number β.So we turn to discuss this dependence next.

V. SENSITIVITY TO β
This section deals with the considerable sensitivity of the model to the scaled contact number β.First, we describe the mathematical aspects, then we remark on the possible disease control implications of this model sensitivity.We note that we already performed a similar study in [2] for the city of Riyadh, and it was found that the system was very sensitive with respect to β.Here, we perform it for the whole country of Saudi Arabia and very similar results are obtained.For the sake of completeness, we do it for the city of Riyadh too.This sensitivity may have considerable policy implications.
We selected three typical examples that predict very different scenarios, with very close values of the contact rate β, and we run the simulations for over seven year, actually 2600 days (since the beginning of MERS in Saudi Arabia).We used the   I), and the two additional values β = 0.1824 and β = 0.1830.The choice was such that the stability numbers were R c < 1, R c = 1 and R c > 1, respectively, but the first and last with values very close to 1.The simulations are depicted in Fig. 7, where the case with β = 0.1824 is depicted in solid curves, β = 0.1824 in dashed curves, and β = 0.1830 in dash-dot curves.The predicted cumulative cases of MERS at the end of the seven years (T) were found to be about 2200, 4200 and 9580, respectively.The cumulative deaths were found to be about 760, 1450 and 3280, respectively.A noticeable differ-ence among the three cases was found, indeed, the numbers more than doubled from the first to the second and from the second to the third case, while the difference between consecutive values of β was just 0.3%.
This clearly indicates that the model is very sensitive to the value of β.As was pointed out in [2], it is very unlikely that this just a mathematical model.This belief is also supported by the description in the literature on the virulent spread of MERS in confined places.
We conducted a similar study of seven years for the city of Riyadh, with results depicted in Fig. 8.We used three scenarios with the values β = 0.1222 (solid lines), β = 0.1231 (dashed lines), and β = 0.1240 (dash-doted lines).The choice was based on the same considerations as for the whole country.The results were about 850, 1740 and 4442 cumulative cases of MERS, respectively, shown in Fig. 8 (T); and about 265, 540, and 1365 cumulative deaths, respectively, shown in Fig. 8 (B).
It is seen that changes of 0.7% lead to the doubling of the cumulative cases of MERS and of the cumulative numbers of deaths.Again, this reinforces the issue of the extreme sensitivity of the model to the scaled contact number.

VI. THE MODEL ERRORS
In this short section we provide a graphic representations of the errors, which are the differences between the data points and the model solution results.These are depicted in Figs. 9 to 11.In Fig. 9 we show the difference between the cumulative reported cases of MERS in Saudi Arabia, for the 865 days used to find the system coefficients using the optimization routine in MATLAB.It represents the errors in the results depicted in Fig. 2 above.Next, Fig. 10 presents the errors for Saudi Arabia in the cumulative numbers of reported cases and the deaths.These are the details presented in Fig. 3 above.Finally, the errors between the data and model simulations for the city of Riyadh in the cumulative numbers of reported cases and the deaths are depicted in Fig. 11.These are taken  from the results in Fig. 5.It seems that the errors do not have any noticeable pattern, with mean about zero, which reinforces our confidence in the model predictions.

VII. CONCLUSIONS
This work deals with the possible trajectories of the MERS disease in Saudi Arabia and in the city of Riyadh.It is a continuation of the study in [2] where the basic model was constructed and simulations of the outbreak of MERS in Riyadh were conducted.
Our aim in this work was two-fold.First, we established the local stability of the DFE and the global stability of the both the DFE and the EE by using only the effective reproduction number or stability control number R c , and this improves the theoretical results in [2], Section 3.3.The analysis of the stability of the model's equilibrium points can be found in Section 3.
The second and more important aim was to simulate and predict the disease outbreak in Saudi Arabia in the very near future, actually, the next two years.We fit the model parameters to a part of the available data, from the first 865 days since the disease was identified, to see how does the model compare with the data for the next 683 days for which data is available, and then used the model to predict the disease outcomes for the next three years, assuming that its trajectory remains the same.It was seen that for both Saudi Arabia and the city of Riyadh, the model predictions for the last 683 days were excellent.Nevertheless, considering the future predictions of the model some caution is in order.Indeed, the baseline simulations for Saudi Arabia, which agreed very well with the data for the 1550 days since the disease was identified, were with the control number R c = 0.99704.However, another parameter fit, which was almost as good, was with R c = 1.004.In the first case there was no endemic equilibrium (EE), while in the second case the EE was found to be asymptotically stable, and these explain why the predictions for the next three years somewhat diverge.In the first case the model predictions for the next three years, until Nov. 4, 2020, there will be about 2200 new cases, the cumulative recovered will be about 1449, and the cumulative deaths will be around 760, Fig. 3.In the second case the model predictions were over 4000 new cases, the cumulative death will be around 2000, and there will be about 2000 recovered.The difference in the predictions is noticeable, although in a country of 32 million population these do not seem to be too divergent.
However, at this stage of the research it is not clear which scenario will play itself in the long run, the one with R c = 0.99704 in which the disease dies (although in 20 years or so) or the one with R c = 1.004, in which the disease is endemic and lingers for a long time.Nevertheless, both predictions seem very reasonable and only time will tell which would be closer to field observations.We stress again that these observations depend crucially on the assumption that the disease continues its current trend.
Similar observations were found for Riyadh, provided in Section 5.
One of the main mathematical features of the model, already pointed out in [2], is its considerable sensitivity to the value of the scaled contact number β.The number measures the probability that one contact between a susceptible individual and a sick one results in infection, and therefore it includes the rate at which people meet each other.Indeed, as was seen in Section 5, Figs.7Biomath 7 (2017), 1802277, http://dx.doi.org/10.11145/j.biomath.2018.02.277 and 8, for Saudi Arabia and for Riyadh, changes of 0.3% and 0.7%, respectively, led to quite different predictions, increases of more than 100%.This sensitivity may have considerable policy implications.In settings where many people congregate and contact is high, the value of may be β higher with possibly severe outbreaks of MERS.

Fig. 1 :
Fig. 1: Compartmental structure and flow chart for the model

Fig. 2 :
Fig. 2: MERS model parameters fit to daily reported cumulative new cases data -red pointsobtained from Saudi Arabian Ministry of Health website during the first 865 days of the disease outbreak.The solid blue line represents the baseline model prediction.The estimated parameters are provided in TableI.

Fig. 3 :
Fig. 3: Saudi baseline simulations of cumulative cases of MERS (T) -green curve; cumulative number of recovered (M) -blue curve; and cumulative number of death (B) -brown curve.The red dots are the field data.The run was for 2555 days (∼ 84 months).

Fig. 4 :
Fig. 4: Saudi simulations with R c = 1.004.Cumulative cases of MERS (T), cumulative number of recovered (M), and cumulative number of death (B).The red dots are the field data.The run was for 2555 days (∼ 84 months).

Fig. 5 :
Fig. 5: Riyadh baseline simulations of cumulative cases of MERS (T) -green curve; cumulative number of recovered (M) -blue curve; and cumulative number of death (B) -brown curve.The red dots are the field data.The run was for 2555 days (∼ 84 months).

Fig. 6 :
Fig. 6: Riyadh simulations with R c = 1.0045.Cumulative cases of MERS (T), cumulative number of recovered (M), and cumulative number of death (B).The red dots are the field data.The run was for 2555 days (∼ 84 months).

Fig. 7 :
Fig. 7: Saudi Arabia.Simulation results of the first case, for the cumulative reported cases of MERS (T) and the cumulative number of deaths (B).Parameter values used are β = 0.1818 -solid curves; β = 0.1824 -dashed; and β = 0.1830dash-dot curves.

Fig. 9 :Fig. 10 :Fig. 11 :
Fig.9: The difference between the data and the model solution in Fig.2

TABLE I :
Model baseline parameters for Saudi Arabia and Riyadh

TABLE II :
Model parameter for Saudi Arabia, the case R c = 1.004.

TABLE III :
Model parameter for Riyadh, the case R c = 1.0045.