Modeling Epidemics as First-order Systems – COVID-19 Example

The semi-logarithmic plot of the cumulative number of cases of epidemics resembles the response of a first-order systems for a step load. This similarity was utilized to develop a first order model that can be used for extracting information about the dynamics of infectious disease epidemics. The developed model was validated using COVID-19 data of China. It was also heuristically fitted to other 13 countries. Obtained results indicated that the model can reliably forecasts the number of infected person, epidemic growth speed towards steady-state condition (process time constant, T), and time to reach steady-state condition (4T). The developed model will help public health authorities in developing more effective control strategies of epidemics.


Introduction
Mathematical models are useful tools that help provide insights into infectious disease dynamics (Caccavo 2020). They can be utilized to get vital information about dynamics of epidemics such as epidemics growth rate, time to reach the peak and the number of infected persons. Such information is of great importance for public health and other concerned authorities. It is useful to develop appropriate control measures and be prepared to deal with consequences of epidemics (Rabbani et al 2020). Due to increase in epidemics occurrence, use of mathematical models has increased tremendously in the last decades. In fact, mathematical modeling has become the main tool for developing effective control strategies for infectious disease epidemics .
Large numbers of mathematical models have been developed to describe the growth patterns of infectious disease epidemics (Lin et al. 2020). Epidemiological models can be divide broadly into statistical and dynamical models. The commonly used dynamical models are either logistic growth models or compartmental models. The first logistic model was proposed by Verhulst in 1834 and modified by Richards in 1959 (Batista, 2020). On the other hand, the first compartmental SIR (susceptible, infected and recovered model) model was developed by the Scottish Kermack and McKendrick in 1927 (Sanglier et al. 2020). Since then, many new versions of both the logistic growth and compartmental models have been developed (Batista 2020; Liu 2020; Ma 2020; Roda 2020). Despite that, however, the original versions of these models are still in use. As a matter of fact, they have been widely used for modeling COVID-19.
During epidemics there is always a need for a fast and effective response (Sanglier et al. 2020). Effective response, however, requires reliable mathematical models that not only fit the profiles of epidemics, but also reliably explain their dynamics. Unfortunately, results of most of the existing models are unreliable because their parameters have no physical meanings and they cannot be accurately estimated from data. Roda et al. 2020 reported that forecast reliabilities of epidemiological models are greatly affected by the identifiability problem of their parameters. Abusam et al. (2020) have shown that parameters of Verhulst and Richards models cannot be identified from data and their predictions are unreliable. According to Louchet 2020, many countries could not reliably predict the progress of the on-going COVID-19 pandemic.
For developing effective control strategies of epidemics, there is always a need for a reliable model that can accurately capture the dynamics of the epidemics (Chowell et al. 2016). In this paper, principles of systems engineering were used to develop a new dynamical model that extracts accurate and physically meaningful information about the dynamics of infectious disease epidemics.

Data
Data used in this paper was obtain from Worldometers website (Worldometers, 2020). Except for China, the data is the daily records of the cumulative confirmed cases, from the 15 February to 30 June 2020. Notice that the data of China used here is for the first wave of COVID-19 epidemic which occurred between 22 January to 2 March 2020.

Model development
The semi-logarithmic plot of cumulative cases of an epidemic grows asymptotically towards a steady state condition (Fig. 1). Data plotted in Fig. 1 is for COVID-19 in China (Worledimeters, 2020) and SARs, which occurred in Singapore in 2003 (Ang, 2007). The two curves shown in Fig.  1 resemble the response of first order systems to step loads. Accordingly, it has been assumed that the logarithm of cumulative cases of an epidemic can be modelled as a first order system subjected to a sudden step load (i.e. increase in number of infected people). Therefore, important information about the dynamics of epidemics (e.g. time constant and time to reach steady-state condition) can accurately be extracted from the data. Notice that conversion of data from arithmetic scale to logarithm scale will just condense it but will not change it. Further, any information obtained in logarithmic scale can easily be converted back to arithmetic scale through calculation of the antilogarithm. Furthermore, all information about time is and will remain in arithmetic scale. First order systems are modeled by a first order differential equations. Denoting for the logarithm of the cumulative cases by C' (pronounced C prime), the response of the system to a step input load u (a sudden change in number of infections) can thus be described mathematically as follows: Where C': Logarithm of the cumulative number of infected cases.
T: Process time constant (day). u: Step input load (sudden change in number of infected persons).
k: process gain (increase in number of cumulative infected persons).
u: process input load (sudden increase in number of infected persons).
Solution of Eqn. 1 is (Bequette, 1998): Where du: Change of input load w.r.t. initial steady state conditions.
Notice that it is not necessary that the initial steady-state be equal to zero (no infected person). Initial steady-state greater than zero means that certain number of infection was sustained before the epidemic appeared.
Assuming the product of k and du equals K, Eqn. 2 can be re-written as: After rearranging and taking the natural logarithm, Eqn. 3 can also be re-written as a straight line equation: Hence for first-order system, the plot of ln[K-C'] is linear (Eqn. 4). To test that, ln[K-C'] of data for Kuwait and Oman was plotted and presented in Fig. 2. From this figure, it is apparent that the plots are not perfectly straight lines, but they can reasonably be approximated by straight lines, except for the first day(s) of the epidemic. Since modeling is a simplification of reality, it is acceptable to assume the plots of ln[K-C'] is a perfect straight lines, when the first day(s) is excluded. That means the initial steady-state condition before the epidemic, C'(0), should be considered an unknown parameter to be estimated from data. This assumption seems to be reasonable since the first record(s) of data is often not reflecting the true number of persons initially infected. With unknown initial steady-state condition, the model will be: The developed model (Eqn. 5) has only three unknown parameters (C'(0), K and T) to be estimated from data. The parameter C'(0) is the number of infected persons just before the epidemic become apparent. The parameter T is the speed at which the epidemic grow (time constant of the epidemic, T) to reach the steady-state condition. From the value of T, time will be taken to reach steady-state condition can also be estimated. Since first order systems achieve more than 98% of their ultimate gain in time equals 4T, therefore they are considered reached steady-state at 4T (Bequette, 1998). Thus, accurate estimation of these physically meaningful parameters will provide public health decision makers with very important information about the dynamics of epidemics.

Model Calibration and Validation
Data of the first wave of COVID-19 epidemic in China was mainly chosen to test the model (Eqn.

5) because it had already reached its end and thus all actual information about it is available.
For calibration, the model was run with data of the first 20 days only. The least square technique, through MATALB function lsqcurvefit, was used to fit model's prediction to the data. Calibration of the model had resulted in estimation of the values of the three unknown parameters: C'(0), T and K.
To validate the model, it was then run with whole 41-days data, using the same parameters' values obtained during the calibration step.
Dynamics of epidemics usually differ from a community to another and thus from a country to another. To demonstrate the ability of the model in identifying the different dynamics of the same epidemics in different countries, the model was also heuristically fitted to COVID-19 data of the following 13 countries: Bahrain, France, Kuwait, Italy, Netherlands, Oman, Qatar, Singapore, Spain, Saudi Arabia, UAE, UK, and USA.

Results and Discussion
Values obtained from model calibration for the three parameters C'(0), T and K, were 5.7296 ± 0.0027, 9.6455 ± 0.0149 and 5.6893 ± 0.0022, respectively. Notice that values of C'(0) and K are in logarithmic scale, whereas that of T is in linear scale. It should also be noted that the standard errors, which were estimated using the Jacobian matrix (Abusam et al 2000 and2001), are much smaller than parameters' values. These results clearly indicate that the three parameters can be identified accurately and reliably from the data.
The anti-log of the value estimated for C'(0) (5.7296) is 308. This means that there were initially 308 persons infected before the step increase in number of infections which resulted in the appearance of the epidemic. Given that the anti-log of the prediction for the first day is 520 (571 from data). The step increase in number of infection is thus 212 persons (520 -308). Therefore, the sudden increase in number of infected persons by 212 had caused the epidemic to appear.
From model calibration, the estimated process time constant (T) is 9.65 days. This means that the epidemic will reach steady-state condition in about 39 days (4 x 9.65 days). Therefore, that is a highly reasonable estimate since the epidemic in China was considered to be over in 41 days. The epidemic in China finished so fast because very strict control measures had been implemented (Lin et al 2020). In short, results of model calibration seemed to be very realistic and accurate.  Results of the heuristic fitting to the data of 13 countries other than China are given in Table 1.
The small values of RMSE, which ranges from 0.0881 to 0.4504, indicate that the model fits very well the dynamics of COVID-19 at the different countries. This can also be seen from Fig. 4 which presents the model fit to the data of Italy as an example. The highest value of RMSE (0.4504) is obtained for Singapore. As can be from Fig. 5, the fit of model with data of Singapore is not as good as that, for example, for Italy (Fig. 4). The relatively large value estimated for the time constant of Singapore (93 days) indicates the epidemic there still has long time (day 376) to go before approaching its end. In spite of the relatively poor fit, however, all the model parameters were also accurately estimated from the data of Singapore as indicated by small standard errors (Table 1).    (Table 1). This demonstrate clearly that the model is able to predict the difference in dynamics of the epidemic under different conditions. Table 1 shows that the values estimated for the initial steady-state conditions (C'(0)), process gain (K) and the process time constant (T) are unique for each country. However, it can be noticed that the initial steady-state conditions (C'(0)) for the European countries (France, Italy, Netherlands and Spain) is relatively smaller (15 -20 days) than that for the other countries, e.g. for Kuwait (109 days). This means that the model predicted that the epidemic in the four mentioned European countries will reach steady-state conditions in only 60 to 80 days, while it that will take 435 days in Kuwait. This seems to be true. The four European countries have already started to ease the control measures and to return normal life, while Kuwait is still implementing some strict regulations.
It becomes clear from above that the developed model can be used reliably for extracting the dynamics of COVID-19 epidemic. However, the same model can also be used to model other epidemics. As shown for SARs in Fig. 1, the semi-logarithmic plots of the epidemiological curves look like a response of a first-order system to a step load. Therefore, the same model can also be used to modeled other epidemics.

Conclusions
A simple first-order system model, with physically meaningful parameters, has been developed for infectious disease epidemics and illustrated using the data of COVID-19 pandemic.
Model calibration has indicated that all the three unknown parameters of the model (state initial condition, process gain and process time constant) can be accurately identified from data.
The developed model has also proved to be reliable at revealing the difference in the dynamics of COVID-19 pandemic in other 13 countries.