Mathematical model and simulations of MERS outbreak : Predictions and implications for control measures

The Middle East Respiratory Syndrome (MERS) has been identified in 2012 and since then outbreaks have been reported in various localities in the Middle East and in other parts of the world. To help predict the possible dynamics of MERS, as well as ways to contain it, this paper develops a mathematical model for the disease. It has a compartmental structure similar to SARS models and is in the form of a coupled system of nonlinear ordinary differential equations (ODEs). The model predictions are fitted to data from the outbreaks in Riyadh (Saudi Arabia) during 2013-2016. The results of model simulations indicate that MERS will eventually be contained in the city. However, the containment time and the severity of the outbreaks depend crucially on the contact coefficients and the isolation rate constant. When randomness is added to the model coefficients, the simulations show that the model is sensitive to the scaled contact rate among people and to the isolation rate. The model is analyzed using stability theory for ODEs and indicates that when using only isolation as the control, the endemic steady state is locally stable and attracting. Numerical simulations with parameters estimated from the city of Riyadh illustrate the analytical results and the model behavior, which may have important implications for the disease containment in the city. Indeed, the model highlights the importance of isolation of infected individuals and may be used to assess other control measures. The model is general and may be used to analyze outbreaks in other parts of the Middle East and other areas. Keywords-Middle East Respiratory Syndrome (MERS); compartmental continuous model; stable endemic steady state; random coefficients; simulations;


I. INTRODUCTION
Middle East Respiratory Syndrome Coronavirus (MERS-CoV) emerged as a new human coronavirus in June 2012.The first MERS patient reported in Jeddah, Saudi Arabia, had pneumonia and renal failure with unknown coronavirus (CoV) that was isolated from his sputum [23].Following identification of the new virus, a subsequent case from Qatar, that was being treated in the United Kingdom, was identified.However, the onset of Citation: Nofe Al-Asuoad, Sadoof Alaswad, Libin Rong, Meir Shillor, Mathematical model and simulations of MERS outbreak: Predictions and implications for control measures, Biomath 5 (2016), 1612141, http://dx.doi.org/10.11145/j.biomath.2016.12.141 the new disease was traced back to an even earlier event.Indeed, in April 2012 a cluster of what were diagnosed then as 'pneumonia' cases had occurred in a hospital in Jordan [2], [9].Currently, this infection appears to be geographically linked to the Middle East, with cases that originated in Jordan, Saudi Arabia, Qatar and United Arab Emirates.Since then, it spread to several other countries: Egypt, Kuwait, Turkey, Oman, Algeria, Bangladesh, Indonesia, France, Germany, Italy, the United Kingdom, South Korea and the United States [21].
The virus is classified as coronavirus which is a single stranded RNA virus.It is similar to Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV), which killed almost 10% of the affected individuals between 2002 and 2003.The mode of transmission to humans is not completely defined yet, but evidence suggests that MERS-CoV can be transmitted to humans via both animals and humans, and it spreads among people who are in close contact.Transmission from infected patients to healthcare personnel has also been observed.Camels and bats are currently known as natural hosts and evidence of cases of camel to human transmission has been reported, see e.g., [19].The incubation period for the disease is 2 to 14 days [14], [21].Confirmed cases include only those with a positive polymerase chain reaction in accordance with the laboratory guidelines for virus genetic material.Probable cases are those that have a link with a confirmed case and a clinically compatible illness but without definitive laboratory confirmation.The majority of patients experienced severe respiratory disease while minority were reported to have non-severe disease and some cases were reported as asymptomatic.Severe disease was defined as admission to an intensive care unit, use of extracorporeal membrane oxygenation, mechanical ventilation and vasopressors.The signs and symptoms include fever, cough, sore throat, chills and shortness of breath.Some patients progress to dyspnea and pneumonia.Complications of the disease include severe pneumonia with acute respiratory distress syndrome, septic shock and multi-organ failure that result in death.The disease has a high level of case fatality.Fatal infections were recorded mainly in immune compromised patients with underlying illness and co-morbidity was frequently recorded [1], [14].Currently, there is no specific treatment for the disease.Its management largely depends on provision of organ support, and prevention of complications.In specific circumstances, additional interventions included empiric use of broad-spectrum antimicrobial agents, antivirals and the addition of antifungal agents to minimize the risk of co-infections with opportunistic pathogens [12], [14].
Mathematical approaches to MERS dynamics were studied recently in several publications.A compartmental model was studied in [4] analyzing epidemiological data on the progression of the MERS-CoV outbreak in April-October 2013 in Saudi Arabia.The model includes community and hospital compartments, and also distinguishes zoonotic (index) transmission from human-tohuman (secondary) transmission.The relative contributions of different cases to the overall level of transmission were evaluated.In [5] statistical methods were used to compare the characteristics of MERS and SARS.Some insights, concerning the nosocomial outbreak in South Korea, from data fitting with the closed-form solution of Richard's model, can be found in [11].The main concern in [8] was to design ways to selectively detect and differentially diagnose only MERS cases among imported cases with upper respiratory symptoms.They used a probabilistic model and a Bayesian approach to help the differential diagnosis based on a known incubation period.Given an imported case of MERS, the numbers of new cases and the transmission were estimated in [17].Another compartmental model, based on ODEs, was used in [22] to study the outbreak of MERS in South Korea during May-June 2015.Finally, a discrete compartmental model was constructed and studied in [18], where the aim was to study spatial heterogeneity.
To study the dynamics of the disease and to assess its spread and the effectiveness of various Biomath 5 (2016), 1612141, http://dx.doi.org/10.11145/j.biomath.2016.12.141 measures to contain or eradicate it, this work develops a compartmental mathematical model for MERS, which is based on models for SARS [10] and MERS [4], [22].We choose to use a continuous model and not a stochastic one, although the data was discontinuous and random, since our model captures the dynamics of the diseases well and avoids some of the mathematical difficulties associated with stochastic differential equations.The model describes the dynamics of five populations consisting of susceptibles, asymptomatics or exposed (individuals who carry the virus and can infect others but have no symptoms), infectious individuals with symptoms, isolated and recovered individuals.The virus dynamics and interactions with animals that may carry it are not included in the model at this initial stage of its development.Mathematically, the existence of a unique solution for the system is straightforward, and here the nonnegativity of the solution is shown.It is found that the model has a disease-free equilibrium and may have an endemic equilibrium.Using the system parameters that were obtained by fitting the model to the epidemiological data, it is seen that the endemic state is locally stable and attracting, which means that the disease may linger for a long period of time if treatment or prevention measures are not applied.
A computer program was written in the package MAPLE for the numerical simulations of the model, and the model coefficients were identified by using an optimization program in MATLAB together with data that were collected in the city of Riyadh [20].Therefore, simulations of the model were done for Riyadh.To study the sensitivity of the model to its parameters, simulations with randomness in some of the coefficients were performed.Among the fitted parameters, the system is very sensitive to the contact parameter.This suggests that isolation or quarantining of infected individuals may be an effective way to control the disease spread until possible treatment and vaccination are developed.
We would like to stress that the system parameters were calibrated from the data collected between Nov. 5, 2013 and March 17, 2016, a period of 864 days.The additional data from March 18, 2016 until Nov. 14, 2016, a period of 242 days, was added in the figures and was found to fit very well the predictions of the model with the original parameters, without any need for parameter modifications.This provides a considerable confidence in the predictive capabilities of the model.
Following the Introduction, the mathematical model is developed in Section II.As noted above, it has compartmental structure resulting in a coupled system of five ODEs.Partial analysis of the model is presented in Section III, where the positivity and boundedness of the solution are established.The stability of the two equilibrium states is described in Section IV.Section V describes the numerical algorithm we developed for the numerical simulations of the model, and presents results of some of the simulations.Additional sensitivity analysis and simulations are shown in the remaining sections.In Section VII we point out the importance of the main contact parameter and its influence on disease dynamics.We conclude the paper in Section VIII, where some unresolved issues and future work can also be found.

II. THE MODEL
In this section, we present a basic model of MERS dynamics that is based on the SARS model in [10].To fix ideas, we use the terminology of a city in which the disease is spreading.We deal with whole populations, and do not take into account their spatial spread.We assume that the populations are large enough to justify the use of continuous description, that is, using ordinary differential equations.If the geographical distribution of the disease within the city is found to be important, models with partial differential equations will be needed, which considerably increase the mathematical complexity of the model and are not warranted at this stage of our understanding of the disease.As was noted in the Introduction, we describe the dynamics of five subpopulations: susceptibles S, exposed or asymptomatics E (below we use exposed, as is customary, referring Biomath 5 (2016), 1612141, http://dx.doi.org/10.11145/j.biomath.2016.12.141 to the individuals who carry the virus and can infect others but have no symptoms), infected I, isolated J, and recovered R individuals.These are functions of time, measured in days.
The susceptible population consists of all the individuals who are alive and do not belong to any of the other groups.The sub-population E(t) denotes the current number of asymptomatic individuals who have been exposed to the virus but have not yet developed clinical symptoms of MERS.The sub-population I(t) denotes individuals with clinical symptoms of MERS.The sub-population J(t) denotes isolated symptomatic individuals, and the sub-population R(t) denotes the individuals who recovered from MERS.It is assumed that those who are exposed become sick before recovering and those who are recovered do not become sick again.These are model assumptions and at this stage we do not know if this is the case, since there is no available data.
The disease behavior that is reflected in the model is as follows.MERS-CoV case fatality rate is very high, 40%, which explains why each suspected case of MERS-CoV is isolated immediately, and when the PCR test is positive for the coronavirus, the patients are isolated in a negative pressure room to limit the spread of the infection.The median incubation period is approximately five days (range 2-14 days) and the median time from illness onset to hospitalization is approximately four days.In critically ill patients, the median time from onset to intensive care unit is approximately five days, and median time from onset to death is approximately 12 days.The recovery is observed to be different from one person to another and depends on the health status and the approximate time to recovery ranges from one to two weeks if there are no complications and extends to a month or more if complications develop.
The flow diagram of the model is depicted in Fig. 1.
We denote by P (t) the number of individuals that are added to the city by birth and from the outside each day, and the natural death rate (in absence of MERS) of the population is µ (1/day).The probability that one contact between a susceptible and an infected results in infection is β (1/day); the probabilities that one contact between a susceptible and an asymptomatic or an isolated results in infection are E β (1/day) and J β (1/day), respectively, where 0 ≤ E , J ≤ 1 are the reduction parameters in the infection rates.
The parameter k (1/day) denotes the rate of development of clinical symptoms in exposed individuals.It is assumed that there is no diseaseinduced death in the exposed class.We denote by γ (1/day) the isolation rate of infectives, which is a control variable since an increase in the isolated may substantially reduce the spread of MERS if done properly.The rate coefficients d 1 and d 2 , with units of 1/day, are the disease-induced death rates of the infected and isolated, respectively, which are in addition to the natural death rate µ.Finally, σ 1 and σ 2 are the recovery rates of the infectives and isolated, respectively.These model assumptions lead to the following mathematical model for the dynamics of MERS.
Model 1: Find five functions (S, E, I, J, R), defined on [0, T ], such that for 0 < t ≤ T , the following hold: The initial conditions are (7) Here, S 0 > 0 is the initial population before the break-out of the disease, and E 0 , I 0 , J 0 and R 0 are nonnegative initial populations.In practice, one typically assumes that S 0 = N (0) > 0, so that at the start of the epidemic there are only susceptibles, and all the other populations vanish.However, for the sake of generality we allow nonnegative initial populations and, moreover, N (0) = S 0 +E 0 +I 0 +J 0 +R 0 so that (1) holds at t = 0.A summary of the definitions of the model parameters is given in Table I.
Equation (2) describes the rate of change of the susceptible population S(t).The second term on the right-hand side describes the rate at which the susceptibles become infected with the virus by being in contact with exposed, infectives and isolated individuals.The last term is the natural mortality term.The rest of the equations have similar structure and similar interpretation.It is seen in ( 2) and ( 6) that the recovered do not get infected again, which is one of the model assumptions.
The period during which an individual is exposed or asymptomatic is 2-14 days.The natural death rate is denoted by µ and using a life expectancy of 80 years, we obtain that µ = 0.000034 per day.In the absence of disease, the total population of the city is N = P/µ = 5 million.
Then, the added susceptibles population per day is P = 170.We note that P includes birth and net population movement between the city and the outside.The rest of the parameters' values used in the simulations are given in Table 2.
Finally, the cumulative deaths induced by the disease were obtained from the expression with the initial value D(0) = 0. We turn to the analysis of Model 1.The existence and uniqueness of a local solution follows from standard considerations for ODEs, since the right-hand side of the system is locally Lipschitz continuous.Next, to show that the model makes sense biologically, we show that the solutions are non-negative.
Proof: Let 0 < t < T and denote by λ the Biomath 5 (2016), 1612141, http://dx.doi.org/10.11145/j.biomath.2016.12.141 force of infection, where (1) holds, i.e., N = S + E + I + J + R. It follows from (2) that Using separation of variables and integration on [0, t], together with the initial condition yields since S(0) > 0. Similarly, for (3), we have which yields It follows from the third and fourth equations that I(t) ≥ I(0)e −(γ+d1+σ1+µ)t ≥ 0, and J(t) ≥ J(0)e −(σ2+d2+µ)t ≥ 0. Finally, Hence, R(t) ≥ R(0)e −µt ≥ 0. Thus, the solution is nonnegative whenever it exists.Next, we show the the boundedness of the solution.It follows from the previous result that the solution is non-negative.By adding the five equations of the model and using (1), we get the rate of change of the total living population N (t), Therefore, Hence, dN dt + µN ≤ P, so that We conclude that N (t) is positive and uniformly bounded from above, for all 0 ≤ t < T .Next, since N (t) > 0 and each one of S, E, I, J and R is non-negative, it follows that they all are uniformly bounded, too.This completes the proof.This result leads to the following corollary about the global existence and uniqueness of the solution.Because of its importance, we state it as a theorem.
Theorem 3: The unique solution (S(t), E(t), I(t), J(t), R(t)) exists for all t ≥ 0.Moreover, the solution remains in the set which is invariant and compact.Since the solution is bounded, the system is uniformly Lipschitz continuous and the global existence follows.

IV. STABILITY OF THE EQUILIBRIUM STATES
It is straightforward to see that the model has two steady states: a disease-free equilibrium and an endemic equilibrium.We next study the stability of these two states.We note that most of the algebraic manipulations were done in MAPLE.

A. Stability of the disease-free equilibrium
The disease-free equilibrium is the state S = S 0 = P/µ and E = I = J = R = 0, obtained by setting the right-hand sides of the equations to zero.Moreover, N = S = S 0 , therefore, the model reduces to µS = P , which yields the relationship S 0 = P µ , and the rest of the variables vanish.We turn to finding the basic reproduction number, R 0 , which controls the disease-free equilibrium and the effective reproduction or control Biomath 5 (2016), 1612141, http://dx.doi.org/10.11145/j.biomath.2016.12.141 number R c , which controls the stability of the endemic equilibrium, see e.g., [10].We note that in what follows the total population is kept constant, N = S 0 .
Let x = (E, I, J, R, S) T , then the model can be written in the form where Next, to find both reproduction numbers, we follow [7] and restrict the system to the infected populations.To compute the Jacobian J a of the reduced system, we evaluate the partial derivatives of F with respect to (E, I, J).Thus, Similarly, evaluating the partial derivatives of V with respect to (E, I, J) gives Next, the inverse matrix of V is given by To calculate R c we need the spectrum, actually the spectral radius, of the matrix where Then, the eigenvalues of F V −1 are given by (12) Hence, the control reproduction number is We note that in various works, see e.g., [10] and references therein, the basic reproduction number R 0 is defined in the absence of any control measures, which in our case means γ = 0 and so it is obtained from R c above when we set γ = 0, and used to study the stability of the disease-free equilibrium.Thus, We remark below on the relationship between R 0 and R c .
Here, we assume that γ = 0 so there are no new infected individuals that become isolated, however, since the result is locally near the origin, some isolated individuals might have been present initially.Moreover, for the sake of uniformity, we retain the 5D system.We have the following stability result.
Proposition 4: The disease-free equilibrium of the model is locally asymptotically stable if R 0 < 1 and is unstable when R 0 > 1.
Proof: The Jacobin matrix of the model is given at the disease-free equilibrium, P 0 = ( P µ , 0, 0, 0, 0), by J a (P 0 ) = This matrix has three negative eigenvalues: −µ (double) and −(σ 2 + d 2 + µ).The remaining two eigenvalues are obtained from the 2 × 2 submatrix The trace of A is given by and it is strictly negative.Next, the determinant is given by Dividing by (k + µ)(d 1 + σ 1 + µ) and using (14) we get, Therefore, since the eigenvalues of A are both negative when det(A) > 0 (recall that tr(A) < 0), we conclude that P 0 is locally asymptotically stable if and only if R 0 < 1.When R 0 > 1, the matrix has a positive real eigenvalue and this means that the disease free equilibrium is unstable.This completes the proof of the Proposition.

B. Existence and stability of the endemic equilibrium
We denote the endemic equilibrium of the model by P * = (S * , E * , I * , J * , R * ) such that P * = (S 0 , 0, 0, 0, 0).To find when it exists, we set the right-hand side of each equation of the model to zero.We solve the resulting system in terms of the equilibrium force of infection at steady state λ * , given by where N * = S * + E * + I * + J * + R * .The system is given by From equations ( 16) and( 17), we obtain We rearrange equations ( 18)-( 20) and write I * , J * and R * in terms of E * .Thus, Substituting ( 21)-( 24) in (15) and simplifying yields, Dividing by D 1 D 2 D 3 and recalling (13), we obtain Thus, The equation has two solutions, λ * 1 = 0 which corresponds to the case when R c = 1, and Biomath 5 (2016), 1612141, http://dx.doi.org/10.11145/j.biomath.2016.12.141 which corresponds to the endemic equilibrium, when positive.We conclude that when R c < 1 the only equilibrium is the disease-free one.The endemic equilibrium exists and is unique when R c ≥ 1.
We summarize this in the following result.Proposition 5: When R c > 1 there exists a unique endemic equilibrium P * = (S * , E * , I * , J * , R * ), which is locally asymptotically stable.
Proof: It remains to show the local stability.Let S = x 1 , E = x 2 , I = x 3 , J = x 4 , R = x 5 (note that the order here is different from the one above), denote x = (x 1 , . . ., x 5 ) and then N = x 1 + x 2 + x 3 + x 4 + x 5 .We rewrite the system (2)-( 6) in the form where the components of f and the system are given as The Jacobian of the system (27)-(31) evaluated at the disease-free equilibrium .
Let β = β * be a bifurcation parameter and R c = 1 be the bifurcation point.Recall that the control reproduction number R c is given in (13).Solving the system for β * when R c = 1, we find Actually, β * = β/R c and R c is a linear function of β.
We note that zero is a simple eigenvalue of J a (P 0 ), hence, we may use the center manifold theory to analyze the system (27)-(31) near β = β * .
The Jacobian J a (P 0 ) has a right-eigenvector associated with the zero eigenvalue at β = β * given by w = [w 1 , w 2 , w 3 , w 4 , w 5 ] T , where in terms of w 3 = η > 0, we have ] of J a (P 0 ) associated with the zero eigenvalue at β = β * is given in terms of v 2 = ξ > 0 by The associated non-zero partial derivatives of f at the disease-free equilibrium are: Biomath 5 (2016), 1612141, http://dx.doi.org/10.11145/j.biomath.2016.12.141 To obtain the result, we use item (iv) in [3,Theorem 4.1].To that end, we let Therefore, we have a < 0 and b > 0 and it follows now from item (iv) in [3, Theorem 4.1] that the unique endemic equilibrium point is locally asymptotically stable when R c > 1 (which also means that β * < β).
The force of infection λ * 2 and I * are functions of R c , and R c is a function of the contact rate β.At R c = 1, there is a transcritical bifurcation and the endemic equilibrium at R c > 1 is locally asymptotically stable.
Remark 6: It was found above that the stability of the disease-free equilibrium is related to R 0 , while that of the endemic equilibrium to R c .An attempt to study the stability of the disease-free equilibrium in terms of R c leads to various additional conditions that do not seem to have much merit.However, using R 0 < 1 as the condition for the asymptotic stability for the disease-free equilibrium and R c > 1 as the condition for the asymptotic stability of the endemic equilibrium leaves a gap in the results.As can be seen below, in the baseline-case in Section V, R 0 = 2.033 while R c = 1.009, so the disease-free equilibrium is unstable and the endemic equilibrium is locally asymptotically stable, and our numerical simulations indicate that it is globally asymptotically stable.However, the mathematical issue is still unresolved.
We next discuss the relationship between the two numbers.It follows from the definitions, in terms of the model parameters, that and Then, Therefore, the relationship between R 0 and R c is determined by ρ * .It is seen in (33) that when d 1 and d 2 , and σ 1 and σ 2 have comparable values, and J << 1, then R 0 > R c .Indeed, as was noted above, in the baseline we have R 0 > R c > 1.
We summarize these observations as follows.
Proposition 7: Let ρ * be given in (33).If and there are three possible cases: Biomath 5 (2016), 1612141, http://dx.doi.org/10.11145/j.biomath.2016.12.141 (i) R c ≤ R 0 < 1.The endemic equilibrium does not exist and the disease-free equilibrium is locally asymptotically stable; (ii) R c < 1 and R 0 > 1.The endemic equilibrium does not exist and the disease-free equilibrium is unstable (further analysis is needed); (iii) 1 < R c ≤ R 0 .The endemic equilibrium is locally stable and attracting and the disease-free equilibrium is unstable.If and there are three possible cases: Both the disease-free and the endemic equilibria are locally asymptotically stable (further analysis is needed); (vi) 1 < R 0 < R c .The same as (iii) above.We conclude that the gaps that make further study necessary are in cases (ii) and (v).

V. NUMERICAL SIMULATIONS
To study the dynamical behavior of the MERS model, a numerical algorithm was developed and implemented in MAPLE and extensive numerical simulations were run, using the parameter values listed in Table II.The simulations were run exclusively for the city of Riyadh in Saudi Arabia since the parameters were obtained by fitting the model to the data available for the outbreaks of the disease in Riyadh from Nov. 5, 2013 to Nov. 14, 2016.The parameters P and µ, which are not associated with MERS, are readily available for the city, while all the other parameters were fitted by using an optimization routine in MATLAB, based on the first 864 days.We would like to point out, as was noted above, that the additional data from the following 242 days were found to fit well into the model without any need to change the previously fitted parameters.We fitted the MERS model ( 2)-( 6) to the data (the first 864 days) using MATLAB's lsqcurvefit function, which is part of the optimization toolbox.Cumulative data of the number of reported cases were obtained from the Saudi Arabian Ministry of Health website [20].We fitted the MERS model to the data using initial values for each parameter.Using the optimization routine, we obtained better estimates of the same parameters from the fit.The values of the parameters are given in Table II.
Using the parameter values listed in Table II, the endemic equilibrium values are: 4 Since all the eigenvalues at P * are negative, it follows that P * is locally asymptotically stable.
Our numerical simulations indicate that the endemic equilibrium is globally asymptotically stable, as is seen in Fig. 2, where four trajectories of the system that start at different initial conditions approach P * .However, the proof of the global stability is not available yet.The results of numerical simulation are shown in Fig. 3 with the initial conditions S(0) = 4999990, E(0) = 0, I(0) = 10, J(0) = 0, R(0) = 0.The figures depict the solution with parameters that provided the best fit to the observed cumulative cases of MERS and cumulative number of deaths reported in [20] (in the first 864 days).The observed data, which include additional 242 days, are plotted in red while the model predictions are the smooth colored cures.We note that there were no reported data of the cumulative number of recovered during the first 182 days, so we set it as zero in the middle figure, which explains why the whole red graph is quite below the blue curve.Assuming that the cumulative number of recovered during the first 182 days was the value taken from the graph made the fit much better, but it was decided to provide the graphs with the missing data.
In Fig. 3, the cumulative infected reported cases of MERS are depicted at the top (T), the cumulative number of recovered in the middle (M) and the cumulative number of deaths on the bottom (B).
It is seen that the model with baseline parameter values predicts that if the epidemic continues on its current trajectory, in another 11 months there will be about 1077 reported cases, the cumulative death count from the disease will be about 375, and there will be about 696 recovered.This may turn out to be a very conservative prediction based on the values of β and γ that were obtained by curve fitting.We discuss the model's sensitivity to these two parameters below, since it is found that small changes in β can lead to large changes in the disease dynamics.

B. Simulations with different values of γ
Since the isolation rate γ is currently the main control parameter, we performed simulations with three different values of γ to find out the effectiveness of the infected individual's isolation on the spread of the disease.This also indicates the model's sensitivity to γ.
Fig. 4 depicts the variation of cumulative reported cases of MERS, cumulative number of recovered, and cumulative number of deaths for different values of the isolation rate γ.It is found that as the isolation rate γ increases, the three populations decrease.Indeed, by differentiating R c with respect to γ, we obtain which means that R c decreases as γ increases.This indicates the impact of isolation in reducing the reproduction number.This is clearly seen in Fig. 4 where the lowest curve corresponds to the highest isolation rate γ = 0.17 and the highest curve to the lowest isolation rate γ = 0.1501.
The dependence of R c on γ is shown in Fig. 5, and once R c < 1 the endemic equilibrium ceases  to exist.Moreover, in the absence of isolation (γ = 0), we have that R c = R 0 .Next, following [15], we define γ c , the critical value of γ, by letting Fig. 4: Model simulations with three values γ = 0.1501, 0.16, 0.17 showing the cumulative reported cases of MERS (T), cumulative number of recovered (M), and cumulative number of death (B).The rest of the parameter values are given in Table II.
We note that in the baseline simulations γ c = 0.1540.
We conclude that when isolation is the only disease control measure, the isolation rate constant γ needs to be well above γ c to make MERS die out.However, in the baseline simulation γ = 0.1501 < γ c , and hence the endemic equilibrium is locally asymptotically stable.

VI. PARAMETER RANDOMNESS AND SENSITIVITY ANALYSIS
In this section we introduce randomness into the model parameters.This allows us to better understand the model's dependence on the parameter values.Moreover, to use the model as a predictive tool, we must find out how these predictions change when the values of the model parameters change.Clearly, small changes in the solution caused by changes in the value of a parameter indicate that there is low sensitivity to the parameter and an approximate value is sufficient, while considerable changes in the solution caused by small changes in a parameter indicate that a more precise parameter value is needed for obtaining reliable predictions.We note that the measurability of the solutions with respect to the random parameters is a consequence of the general measurability results in [13].
We introduce randomness into three parameters: βthe effective contact rate, E -the reduction in the transmission rate from exposed to infected with clinical symptoms, and γthe isolation rate of those with clinical symptoms.The rest of the parameters were kept at baseline values.
For each parameter, we let the probability space be (Ω, F, P ), where Ω is the sample space, F is the Borel σ-algebra, and P is the probability function.The latter was chosen in all cases as the uniform probability.We run 100 simulations, each over four years, using the relevant coefficient that was varying randomly in Ω.At each computational time step we choose the highest and the lowest values of the 100 solutions and constructed an envelope that contains all the solutions.These envelopes are depicted in the figures.However, we note that these envelopes themselves are not solutions.In each case we tried to adjust the sample space ω so that the available data (from the first 864 days) falls within the envelopes.
The data were collected for 3 years for the infected and deaths, and two years and 6 months for the recovered (the recovered data were not available in the first 182 days).These simulations provide the predictions for the next 11 months, till the end of October 2017.The last data point was from Nov. 14, 2016.As was mentioned above, the red plotted points are the data while the colored curves are the model predictions.

A. Sensitivity with respect to β
We begin with the system sensitivity to the variation in the effective contact rate β.In the baseline setting β = 0.1220.We let β, the random variable for the contact rate, be given by β = β + ω, where β = 0.1220, and ω was chosen randomly from the sample space We used a random β in each one of the 100 simulation runs, and each simulation was for a period of four years.The choice of Ω β was based on the fact that the values 0.09 ≤ β ≤ 0.1241 led to envelopes that contained the known data.
The results with random β are shown in Fig. 6.The top curve depicts the maximum value and the bottom curve shows the minimum value of the solutions at each day.It is seen that the minimum and the maximum values agree well with the data of the cumulative reported cases of MERS, cumulative number of recovered and cumulative number of deaths.The variations in the cumulative reported cases of MERS are depicted in Fig. 6 (top) where the upper and lower envelopes reach at the end of the fourth year the maximum and minimum values of 1909 and 29, respectively.The cumulative number of recovered (middle) shows maximum and minimum values of 1223 and 19, respectively.The variations in the cumulative number of deaths (bottom) show maximum and minimum values of 659 and 10, respectively.However, the real interest in this simulation is in the predictions of the upper envelopes, which depict the maximal values of the cumulative variables.We conclude that the system is sensitive to the contact rate coefficient β.Below, because of its importance, we discuss the effects of a wider range for β.

B. Sensitivity with respect to E
We turn to the system sensitivity to the transmission reduction factor E .In the baseline setting E = 0.2996.We let E = E + ω, where E = 0.2996 and ω is chosen randomly from the sample space Therefore, E ∈ [0.01, 0.318].Again, the choice of Ω E was such that the resulting simulations include all of the data points (for the first 864 days).
We used a random E in each year of the 100 runs (each one was for four years!).The simulations are shown in Fig. 7.The variations of the cumulative reported cases of MERS are depicted in Fig. 7 (top) where the upper and lower envelopes reach maximum and minimum values of about 2182 and 30, respectively.The cumulative numbers of recovered (middle) reach maximum and minimum values of about 1396 and 19, respectively, and the cumulative numbers of deaths (bottom) reach maximum and minimum values of about 752 and 10, respectively.We find that the model is sensitive to the variations in the transmission reduction factor E .Since the lower Biomath 5 (2016), 1612141, http://dx.doi.org/10.11145/j.biomath.2016.12.141 limits have already been exceeded, we use them to depict the sensitivity to the parameter.Fig. 7: Variations in the transmission reduction factor E lead to variations in the cumulative reported cases of MERS (T), cumulative number of recovered (M), and cumulative number of death (B), as ω varies randomly in Ω E , so that E varies in [0.01, 0.318].

C. Sensitivity with respect to γ
We study next the system sensitivity to the isolation rate γ, which was assumed to be constant in previous simulations, at the baseline value.
We let γ be given by γ = γ + ω, where γ = 0.1501 (the baseline value) and ω is chosen randomly from the sample space Therefore, the isolation rate constant γ varies in the range 0.144 ≤ γ ≤ 0.167.
We used random γ in each one of the 100 runs.The results of the numerical simulations with random isolation rate γ are shown in Fig. 8.It is seen that the minimum and the maximum values agree well with the data (the full 1107 days) for the cumulative reported cases of MERS, cumulative number of recovered, and cumulative number of death in the time period for which the data is available.The variations of the cumulative reported cases of MERS are depicted in Fig. 8 (top) where the envelopes reach maximum values of about 2537 and 238, respectively.The cumulative numbers of recovered (middle) reach maximum values of about 1620 and 155, respectively, and the cumulative numbers of deaths (bottom) reach maximum values of about 876 and 83, respectively.

D. Combined sensitivity with respect to β, E and γ
Finally, for the sake of completeness, we present the system sensitivity to the combined randomness in β, E , and γ.It was assumed that each year the variables (β, E , γ) were given by (β, E , γ) = (β, E , γ) + (ω β , ω E , ω γ ), where the variations were the same as the individual cases above.11, respectively.Again, we find that the system is quite sensitive to the combined random variation of β, E and γ within the indicated ranges.
We note that the results of the combined variations are very similar to those for separate variations.

E. Comparison of the sensitivity analyses
A comparison of the results above for β, E and γ indicates that within the chosen ranges, the system is sensitive to each one of these parameters, and with comparable sensitivity to all three combined.This reinforces the importance of the contact number and the isolation rate in the model predictions.However, as we discuss in the next section, when we slightly increase the randomness range for the contact rate constant β, the sensitivity result changes considerably.We conclude that the model predicts, assuming that the same trend continues and there are no changes in the virus transmission or new countermeasures, that at the end of the next 11 months (i.e., by end of October 2017) the cumulative cases of reported, recovered and deaths will be below 2320, 1480 and 800, respectively.

VII. ON THE SENSITIVITY TO β
The contact rate coefficient β plays a significant role in the spread of MERS disease, as is seen in most infectious diseases.The simulations presented in this work are based on data collected in the city of Riyadh in Saudi Arabia, and the baseline value of β reflects it.We used the value β = 0.1220 that was obtained from the best fit to the data pertaining to Riyadh.In the sensitivity analysis we found that the range 0.09 ≤ β ≤ 0.1241 well contained the observed data.
However, during mass gatherings, especially where the human density is large, the situation may be very different.In such events the MERS-CoV virus can be easily transmitted through humanto-human contacts, leading to a higher value of β.So we studied the model predictions with larger values of the contact rate.We found, as we depict below, that in such cases the model predicts substantially higher numbers of infected and deaths.An additional important issue that may compound the situation is when such a mass gathering includes individuals from different locations and possibly from different countries.When those who are exposed travel back to their places of origin, they may spread the disease to other places, see e.g., [6].To assess and illustrate the importance of the contact rate, in additional to the sensitivity results above, we run numerical simulations with different values of β.We found that varying the effective contact rate β, while keeping all the other parameters fixed at the baseline values, has a considerable effect on the population.In Fig. 10, the variation of cumulative reported cases of MERS and deaths are shown for two different values of β.A small increase above the baseline value, leads to a large increase in both the symptomatic population and the cumulative number of deaths.This emphasizes the fact that in the absence of vaccination or effective drugs or treatment, the spread of the disease may be largely controlled by decreasing the contact rate.Fig. 10 depicts two very different scenarios: using the baseline value β = 0.12200 results in about 1530 cumulative reported cases of MERS over 5 years; while the slightly larger value β = 0.12887 results in about 298, 178 newly infected cases over the same time period.Similarly, the baseline value results in about 532 cumulative death cases over 5 years, while β = 0.12887 results in about 100, 000 cumulative deaths over the same time period.It is clear that there is a qualitative difference between the two cases, while the difference between the values of β is just 7%.The question whether this sensitivity is just a mathematical artifact of the model or a real property of the disease is an open question that will be resolved with additional data.However, one must be very careful in treating each outbreak of the disease in crowded events by immediate isolation and reduction of human contacts.

VIII. CONCLUSION AND OPEN ISSUES
This work presents a compartmental model for the outbreak of MERS.It aims at the study of the disease dynamics and assessment of various disease containment measures, since currently there is no treatment or vaccination for the disease.
The analysis of the model shows that it is wellposed and possesses a disease-free state and may posses an endemic equilibrium state.The nonnegativity of the solution is shown.Using the baseline values of the model parameters that were fitted from data for the city of Riyadh, the endemic equilibrium state is found to be locally asymptotically stable, which means that the disease may be contained, but not eradicated.To eradicate the disease, as long as treatment or vaccination are not available, one must increase the isolation rate γ and decrease the contact number β so that the disease-free equilibrium becomes asymptotically stable.We note that in such a case the endemic equilibrium does not exist and the disease will eventually die out.
A computer code was written in the package MAPLE and simulations of the disease dynamics were performed.Based on the data from Riyadh gathered over a period of 36 months (during the years 2013-2016) the baseline model coefficients were identified.Actually, the coefficients were identified from the data of the first 864 days of the epidemic.The data for the additional 242 days was found to fit well into the model thus fitted.The predictions for the next 11 months were obtained, based on this data fitting.It is seen that the disease will affect an increasing number of individuals, predicting a cumulative number (over 4 years, till about October 2017) of about 1077 reported cases, 696 recovered and 375 deaths.However, sensitivity analysis indicates that the system is quite sensitive to the contact number β, the transmission reduction factor E and the isolation parameter γ.Indeed, the system sensitivity to β and E is comparable, while the sensitivity to γ is more pronounced.When the combined sensitivity to the three parameters was conducted, it was found that the combined sensitivity was comparable to that of γ, reinforcing the importance of the isolation rate in the model predictions.
When slightly increasing (by 7%) the value of the contact rate constant β, the result changes considerably and the disease becomes more virulent.We conclude from the sensitivity analysis that under the assumption that the same trend continues and there are no changes in the virus transmission or new countermeasures such as vaccination, that by the end of the next 11 months (i.e., by end of October 2017) the maximal cumulative cases of MERS will be below 1909, the recovered below 1223, and the number of deaths below 659.These results may be of help in policy decisions by the health authorities.However, there must be considerable attention paid to the disease spread because small changes in the rates may yield large changes in the numbers of affected individuals.
The model is quite general and flexible.If the assumptions or the rate constants change because of new data, it is easy to modify the model to include the new developments.Moreover, the virus dynamics and interactions with animals that may carry it are not included in the model at this initial Biomath 5 (2016), 1612141, http://dx.doi.org/10.11145/j.biomath.2016.12.141 stage of its development.It may be of interest to investigate if addition of these topics, which will make the model considerably more complex, improves its predictive power.
An interesting issue that was not studied here, but needs to be added is the dependence of the parameters, especially the contact number β, on the state of the immune system of the susceptible person.There may be very different values associated with a person with healthy immune system and one with a compromised one.However, there is no current data that would allow us to address this issue, and the valu of β was obtained from the currently available data.
We note that the question of the global asymptotical stability of the endemic equilibrium when R c > 1 is unresolved, yet, although our numerical experiments indicate that in the baseline case it is, indeed globally, asymptotically stable (Fig. 2).
We plan to use the model to study the disease dynamics in other localities in the Middle East, as there is considerable growth in the disease and the available data.The relationship between R c and R 0 warrants further mathematical analysis.Also, we will employ statistical methods to analyze some of the data, to make the model easier to use by the authorities and policy makers.

Fig. 2 :
Fig. 2: Stability of the endemic equilibrium.Four trajectories in the S−E plane (upper) and the S−I plane (lower) ending in the endemic equilibrium.

Fig. 3 :
Fig. 3: Baseline model simulations of cumulative reported cases of MERS (T -green curve); cumulative number of recovered (M -blue curve); and cumulative number of death (B -brown curve).The red curves are the observed data for 1107 days.

Fig. 8 :
Fig.8: Variations with respect to γ lead to variations in the cumulative reported cases of MERS (T), cumulative number of recovered (M), and cumulative number of death (B) as ω varies randomly in Ω γ so that γ varies in the interval [0.144, 0.167].

Fig. 10 :
Fig. 10: Simulation results for the cumulative reported cases of MERS (T) and the cumulative number of deaths (B).Parameter values used are β = 0.12200 -dashed lower curves; and β = 0.12887 -solid curves.

TABLE I :
Symbols and description of parameters used in the model P recruitment rate of susceptible individuals β effective contact rate E reduction factor in transmission rate by exposed per day J reduction factor in transmission rate by isolated per day k rate of development of clinical symptoms in exposed population µ natural death rate

TABLE II :
Model parameter values.