Modified SEIR Model

A previously developed model^{Reference He, Dushoff and Day30} for an influenza outbreak is adopted for the purposes of this study. The starting point of the proposed model is the SEIR model. The functions S(t), E(t), I(t), R(t) describe the number of susceptible, exposed, infectious and removed (dead or recovered) people at time t. The total population N = S(t) + E(t) + I(t) + R(t) is assumed to be constant. However, an additional extra class is introduced. P(t) expresses the public perception of risk depending on the number of reported confirmed cases that are infectious, in severe/critical condition or dead. Furthermore, the demography³¹ of Turkey and the course of the disease for different age groups is considered (Bock W, Adamik B, Bawiec M, et al. Mitigation and herd immunity strategy for COVID-19 is likely to fail. 2020. https://doi.org/10.1101/2020.03.25.20043109). These data are summarized in Table 1.

TABLE 1

Demography of Turkey on December 31, 2019, and Severity of Symptoms

Several scenarios are presented in this study. The 4 scenarios serve to introduce the model and show the influence of the various parameters considered. First, results for a constant transmission rate are presented. Second, the influence of the average temperature and then the influence of human reactions in the population in combination with the temperature is presented. An exemplary government response is also considered.

The modeling starts at t = 1 day corresponding to March 11, 2020, which is the date of the first confirmed COVID-19 case in Turkey and ends at August 31, 2020, for demonstration purposes in the first 3 scenarios. The transmission rate in the basic form of SEIR models is typically taken as a constant. However, β(t) in Equation (1) is taken as transmission rate in this study. Temperature T(t) effects described by ξ ≥ 0, which is the strength of the response of transmission rate to temperature variations, government actions like school closures encoded in α(t) and individual reactions given by the perception of risk P(t), which decreases when less people die and increases when people die more often due to COVID-19, as well as the response strength κ ≥ 0 are included in Equation (1)

(1) $${{\it \beta} \left( t \right) = {{\it\beta} _0}{e^{ - \xi T\left( t \right)}}\left( {1 - \alpha \left( t \right)} \right){\left( {1 - \frac{{{P\left( t \right)}}\over {N}}} \right)^\kappa }\hbox{.}$$

The codomain of α(t) is [0, 1] while it is [0, N] for P(t). Thus, increasing government actions α(t) in the second factor in Equation (2) or increasing risk perception P(t) in the final factor in Equation (2) lead to a decreasing transmission rate. The proposed compartmental model is given in Equation (2) by

(2) $$\matrix{ {{d \over {dt}}S(t)} \hfill & = \hfill & {{{\beta (t)S(t)I(t)} \over N},} \hfill \cr {{d \over {dt}}E(t)} \hfill & = \hfill & {{{\beta (t)S(t)I(t)} \over N} - \sigma E(t),} \hfill \cr {{d \over {dt}}I(t)} \hfill & = \hfill & {\sigma E(t) - \gamma I(t),} \hfill \cr {{d \over {dt}}R(t)} \hfill & = \hfill & {\gamma I(t),} \hfill \cr {{d \over {dt}}P(t)} \hfill & = \hfill & {d\gamma I(t) - \lambda P(t),} \hfill \cr } $$

where σ⁻¹ is the mean latent period, γ⁻¹ the mean infectious period, d the proportion of severe and critical cases, λ⁻¹ the mean duration of public reaction and β₀ the original transmission rate. In addition to the temperature effect ξ, the humidity should also play a role as seasonal reality. However, to the best of my knowledge, there is no peer-reviewed quantitative information about this influence on August 04, 2020. Note that the information ξ_h on the strength of the response of transmission rate to humidity variations would lead to an additional factor $$e^{ - {\xi _h}h\left( t \right)}$$ in Equation (1). In the framework of this study, the temperature is strongly correlated to the humidity. Hence, climatic effects are modeled with the temperature. The average temperature distribution of the months between March and June in Turkey are taken for the year 2020. The average temperatures of the months between July and December are taken for the year 2019.^{Reference Björk and Strömberg32} These data are summarized in Table 2. In this way, an attempt is made to make the most accurate analysis possible for the months in 2020. All the accessible data are given in the literature.^{33-Reference Neher, Dyrdak and Druelle37} The modeling parameters are summarized in Table 3.

TABLE 2

Average Temperatures of Certain Months in 2019 and 2020

TABLE 3

Parameters of the Model in Equation (2)

Note that the response strength κ has been chosen as highest comparable value found in a previous publication, where the Spanish flu was analyzed,^{Reference He, Dushoff and Day30} which means that the population is exceptionally careful and very well informed, while α₀ is taken as low as possible due to the harsh social economic results of any action the government takes. Moreover, values zero for ξ, κ, and α₀ are unrealistic and serve demonstrations such that it is possible to see the outcome of no reaction of the population as well as the government. The individual response of the population in Equation (1) is increasing with increasing number of infectious people. However, the response of the government α(t) is taken as a varying function due to published evidence.^{Reference Zhao, Lin and Ran38} It is assumed that the government can take action for 2 mo continuously. The measures are then relaxed for 1 mo and used again for 2 mo. In total, the response of the government is given by

(3) $$\matrix{ {\alpha (t)} \hfill & = \hfill & {{\alpha _0}(H(t - 22\,{\rm{days}})H(82\,{\rm{days}} - t)} \hfill \cr {} \hfill & {} \hfill & { + (H(t - 113\,{\rm{days}})H(174\,{\rm{days}} - t))'} \hfill \cr } $$

where H(t) is the Heaviside step function. Day 22 corresponds to April 01, 2020, 82 to May 31, 2020, 113 to July 01, 2020, and 174 to August 31, 2020.

Equation (2) is rewritten in Equation (4) as

(4) $$\matrix{ {\underbrace {{d \over {dt}}\left( {\matrix{ {S(t)} \hfill \cr {E(t)} \hfill \cr {I(t)} \hfill \cr {R(t)} \hfill \cr {P(t)} \hfill \cr } } \right)}_{\left. { = |y(t)} \right\rangle } = \underbrace {\left( {\matrix{ 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill & 0 \hfill \cr 0 \hfill & { - \sigma } \hfill & 0 \hfill & 0 \hfill & 0 \hfill \cr 0 \hfill & \sigma \hfill & { - \gamma } \hfill & 0 \hfill & 0 \hfill \cr 0 \hfill & 0 \hfill & \gamma \hfill & 0 \hfill & 0 \hfill \cr 0 \hfill & 0 \hfill & {d\gamma } \hfill & 0 \hfill & { - \lambda } \hfill \cr } } \right)}_{ = A}\left( {\matrix{ {S(t)} \hfill \cr {E(t)} \hfill \cr {I(t)} \hfill \cr {R(t)} \hfill \cr {P(t)} \hfill \cr } } \right)} \cr { + {{\beta (t)S(t)I(t)} \over N}\underbrace {\left( {\matrix{ { - 1} \hfill \cr 1 \hfill \cr 0 \hfill \cr 0 \hfill \cr 0 \hfill \cr } } \right)}_{\left. { = |{e_2}} \right\rangle \left. {|{e_1}} \right\rangle }} \cr } $$

Equation (4) is solved by use of the classic fourth-order Runge-Kutta method.

Review of the Runge-Kutta Method

Runge-Kutta methods are widely used in different fields like mathematical physics and engineering for solving ordinary differential equations numerically. For this purpose, the first-order initial-value problem in Equation (5) is considered

(5) $$\frac{{{\rm{d}}}\over {{{\rm{d}}t}}}\left| {y\left( t \right)}\rangle \right. = \left| {f\left( {t\comma\left| {y\left( t \right)} \right.} \rangle\right)}\rangle \right.\comma\left| {y\left( {{t_0}} \right)} \rangle\right. = \left| {{y_0}}\rangle \right.\comma\ t \in \left[ {{t_0}\comma\ {t_{{\rm{max}}}}} \right]\hbox{.}$$

Dividing the interval [t ₀, t _max] into M subinterbals [t _m, t _m+1] with m = 0,1, …, M − 1, t_M = t _max and h = t _m+1 − t_m, integrating Equation (5) over each subinterval, using the mean value theorem, defining $$x \in \left[ {{t_m}\comma{t_{m + 1}}} \right]$$ and approximating $$\left| {f\left( {x\comma\left| {y\left( x \right)} \right.}\rangle \right)}\rangle \right.$$ by a linear combination of the vectors $$\left| {f\left( {{x_1}\comma\left| {y\left( {{x_1}} \right)}\rangle \right.} \right)}\rangle \right.\comma \cdots \comma\left| {f\left( {{x_s}\comma\left| {y\left( {{x_s}} \right)}\rangle \right.} \right)}\rangle \right.$$ leads to Equation (6)

(6) $$\left| {y\left( {{t_{m + 1}}} \right)}\rangle \right. = \left| {y\left( {{t_m}} \right)}\rangle \right. + h\mathop \sum \limits_{j = 1}^s {a_j}\left| {f\left( {{x_j}\comma\left| {y\left( {{x_j}} \right)}\rangle \right.} \right)}\rangle \right.\hbox{.}$$

By choosing different values for the parameters s, a_j and x_j different Runge-Kutta formulations can be obtained. The most widely used form, the classic fourth-order Runge-Kutta method, is obtained for s = 4 in Equation (7)

(7) $$\matrix{ {\left. {|y({t_{m + 1}})} \right\rangle } \hfill & = \hfill & {\left. {y({t_m})} \right\rangle + {h \over 6}(\left. {{k_1}} \right\rangle \left. {2|{k_2}} \right\rangle \left. { + 2|{k_3}} \right\rangle + \left. {2|{k_4}} \right\rangle )o({h^5}),} \hfill \cr {\left. {{k_1}} \right\rangle } \hfill & = \hfill & {\left. {|f({t_m},|y({t_m}))} \right\rangle } \hfill \cr {\left. {{k_2}} \right\rangle } \hfill & = \hfill & {\left. {|F({t_m} + {h \over 2},|y\left. {({t_m})} \right\rangle + {h \over 2}\left. {|{k_1}} \right\rangle )} \right\rangle } \hfill \cr {\left. {{k_3}} \right\rangle } \hfill & = \hfill & {\left. {|F({t_m} + {h \over 2},|y\left. {({t_m})} \right\rangle + {h \over 2}\left. {|{k_2}} \right\rangle )} \right\rangle } \hfill \cr {\left. {{k_4}} \right\rangle } \hfill & = \hfill & {\left. {|F({t_m} + h,|y\left. {({t_m})} \right\rangle + h\left. {|k3} \right\rangle )} \right\rangle } \hfill \cr } $$

Hence, $$\left| {y\left( t \right)} \rangle\right.$$ is determined numerically by iteration.

Application to the Modified SEIR Model

The initial-value problem of our modified SEIR model in Equation (4) can be written as Equation (8)

(8) $${{\rm{d}} \over {{\rm{d}}t}}\left| {y\left( t \right)}\rangle \right. = {\bf{A}}\left. {\left| {y\left( t \right)} \right.} \right\rangle + {{\beta \left( t \right)S\left( t \right)I\left( t \right)} \over N}\left( {\left. {\left| {{e_2}} \right.} \right\rangle - \left. {\left| {{e_1}} \right.} \right\rangle } \right)\comma \hskip-64{|} {{y\left( t \right)} \rangle} = \left( {\matrix{ {83154997} \cr 0 \cr 1 \cr 0 \cr 0 \cr } } \right)\hbox{.}$$

Thus, it is of the same form as Equation (5) and the classic fourth-order Runge-Kutta method can be used. Note that the step size can be chosen freely as time steps, such as days or hours, while smaller steps increase the accuracy of the numerical calculation. However, the computation time also increases with decreasing step sizes.