An epidemic-economic model for COVID-19

1.   Introduction

The coronavirus disease 2019 (COVID-19) led to high morbidity and mortality rates throughout the world and created unprecedented challenges in public health and the global economy. Two years have passed since the World Health Organization (WHO) officially declared COVID-19 as a global pandemic in March 2020, but many nations are still struggling with the spread of the infection. The emergence of several new SARS-CoV-2 variants adds further uncertainty to the progression of the pandemic. As of March 2022, more than 470 million cases and 6 million deaths were reported in 225 countries and territories ^[1].

During the first few months of the pandemic, the rapidity of how the SARS-CoV-2 spread across the world caught most nations unawares and poorly prepared to contain the disease outbreaks and treat the most severe cases ^[2,3]. The situation was compounded by the lack of widespread testing, absence of COVID-19 vaccines, and unavailability of effective treatment modalities ^[2]. Furthermore, since the prevalence of COVID-19 among those with mild or no symptoms can be considerably high, tracking the transmission of the disease had proved difficult ^[4]. Nevertheless, fundamental public health mitigation strategies were implemented which included, but were not limited to, personal and environmental hygienic practices, isolation and quarantine of known cases, contact tracing, social distancing, closure of non-essential businesses and services, and shelter-in-place orders ^[5,6]. Such mitigation strategies appeared effective in slowing the transmission of the virus and reducing the number of cases in communities ^[7,8].

The spread of COVID-19 also caused severe downturn and huge uncertainty for the economy throughout the world ^{[9,10,11,12,13]}. In particular, during the "Stay-at-Home" period (from March to May in 2020) in the US, non-essential businesses were shut down and the US experienced a sharp increase in unemployment as of the end of March 2020, representing a 9.5% unemployment rate ^[14]. This increase, up from 3.5% in February 2020, was just one indication on the breadth of the economic upheaval caused by COVID-19. Meanwhile, the US Department of Commerce reported that total retail sales in March 2020 plunged by 8.7 percent from the previous month, the largest decline in the nearly three decades the government had tracked the data ^[15]. The impact on the global economic and financial system were estimated as being at least a 3% reduction in gross domestic product (GDP), with global financial markets experiencing dramatic instability as seen in the stock markets, assets, and risk markets ^[16].

A large body of experimental, clinical and theoretical studies have been generated to better understand COVID-19 and its consequences, and explore more effective intervention strategies (see reviews ^{[12,17,18,19,20]} and references therein). Particularly, many mathematical and statistical models have been published to study the transmission and spread of COVID-19 ^{[21,22,23,24,25,26,27,28,29,30,31,32,33]}, as well as the effects of COVID-19 vaccination ^{[34,35,36,37,38]}. Meanwhile, a number of computational and quantitative studies have considered the impact of COVID-19 on the economic development. For example, de la Fuente-Mella et al. statistically evaluated the effects of the COVID-19 pandemic on the economy of several countries ^[39]. Chen et al. proposed a network-based epidemic-economic model to estimate the direct impact of labor supply shock to each sector arising from morbidity, mortality, and lockdown ^[40]. Altig et al. analyzed several economic uncertainty indicators for the US and UK before and during the COVID-19 pandemic ^[41]. Eichenbaum et al. combined a canonical epidemiology model with macroeconomic models to study the interaction between economic decisions and epidemics ^[9]. Jena et al. developed a multilayer artificial neural network model to forecast the impact of COVID-19 on the GDP of eight countries including the US ^[42]. Xiang et al. combined economic theory with an epidemiological model to explore the long-term impact of the pandemic on economic growth and the effects of different policy responses ^[43]. Dimarco et al. proposed a wealth transfer model to compute the Gini-index in the presence of a pandemic, with an observation of the emergence of economic inequalities ^[44]. In addition, Auld and Toxvaerd estimated behavioral responses to the global rollout of COVID-19 vaccines and found that countries with earlier and higher vaccination coverage strongly tended to be richer ^[45].

In spite of these studies, a fundamental question remains to be answered: How do the COVID-19 transmission and spread, the disease prevention and intervention, and the economic development impact each other within an interconnected triad? The persistence of the pandemic and the wide-ranging slowdown of the global economy at present, long after the onset of the pandemic, underscore the gap between the complex epidemic-economic mechanisms associated with COVID-19 and our current understanding of these phenomena.

In this paper, we seek to partially address this fundamental question by formulating a new mathematical model to quantify the economic impact of COVID-19 and investigate the interplay between the epidemic progression, the economic growth, and the disease management. We consider both the symptomatic and asymptomatic infections in this model, and incorporate the effectiveness of disease control into the respective transmission rates. Meanwhile, factors associated with economics and epidemiology are both included into the model, where the level of pandemic mitigation is negatively correlated to that of economic development. Additionally, the progression of the pandemic and the evolution of the susceptible, infectious and recovered population groups directly impact the mitigation and economic development levels.

To demonstrate a real-world application of this modeling framework, we study the epidemic and economic situations in the US state of Tennessee in 2020. We utilize the unemployment rate, defined by the number of people who are unemployed as a percentage of the labor force, as the key economic indicator in this work. It is known that a high unemployment rate increases economic inequality, generates redistributive pressures, drives people to poverty, and negatively affects long-run economic growth. Our aim is to study how the unemployment rate would intertwine with the mitigation efforts under the impact of COVID-19. We fit our model to the COVID-19 cases and the unemployment rates reported in Tennessee ^[46,47], and conduct detailed numerical simulation to examine the relationship between the epidemic spread, the disease management, and the economic development.

The remainder of this manuscript is organized as follows. We present the model formulation in Section 2, with detailed mathematical analysis provided in Appendices A and B. We then conduct data fitting and numerical simulation using the epidemic and economic data from the US state of Tennessee in Section 3. We conclude the paper with some discussion in Section 4.

2.   Model formulation

We divide the human population into four compartments: the susceptible ($S$), the exposed ($E$), the symptomatic infectious ($I$), and the recovered ($R$). A distinctive feature of COVID-19 is that asymptomatic and pre-symptomatic infection is common ^[18,48], so that infected individuals could be contagious even during the incubation period. We thus assume that the exposed individuals are capable of transmitting the disease; basically, they are regarded as pre-symptomatic infectious individuals in this study. To incorporate the impact of mitigation policies on economic development, we introduce two additional variables: the mitigation level/effectiveness, denoted by $M$, and the economic development level, denoted by $C$. We normalize the range of $M$ such that $0 \leq M \leq 1$, with $M = 1$ representing the situation with the maximum disease control and $M = 0$ the situation with no disease control at all. Similarly, we normalize $C$ to the range between $0$ and $1$ such that $C = 1$ indicates the maximum level of economic development and $C = 0$ indicates the worst scenario of economic development. We assume that the transmission rates are modulated by the disease control, that the mitigation level is stimulated by the disease prevalence, and that the economic development level depends on the available labor supply. We additionally assume that the disease mitigation level and economic development level are negatively correlated to each other.

The following differential equations represent our epidemic-economic model for COVID-19 transmission dynamics, with descriptions for all the parameters provided in Table 1.

The last two equations in system (2.1) reflect the negative correlation between the mitigation level and the economic development level; i.e., the increase of one variable would tend to decrease the other. This assumption has been motivated by empirical observations of the pandemic management and economic growth in the US. Figure 1 (left) shows the quarterly percentage change of the GDP in the US from 2018 to 2020 ^[49]. In particular, under the impact of COVID-19, the GDP in the first quarter of 2020 experienced a reduction of about 5% from the preceding quarter, and the second quarter GDP plunged about 32%. Such reductions were largely due to the closure of businesses, implementation of stay-at-home orders, and other pandemic mitigation strategies that started in March 2020 and extended to May/June 2020. Afterwards, the third quarter saw a large increase of GDP (about 33%) from the previous quarter, in parallel with the re-opening of businesses, removal of the stay-at-home requirement, and other relaxed control measures. The fourth quarter in 2020 continued the trend of GDP growth with a moderate 4% increase. Figure 1 (right) displays the quantification of the business closure and stay-at-home requirement as two representative mitigation policies, where a higher number indicates a stronger effort, for the US in 2020 ^[50]. We can observe that both policies had the highest strengths in April, which then quickly decreased afterwards, reaching a low point in June/July. For the last two quarters in 2020, the two policy curves stayed at approximately the same levels as that in June, though occasional slight increases (e.g., July–August, November–December) can be noticed. These figures convey the message that a focus on pandemic management would slow down economic development, whereas an emphasis on economic growth would necessitate policy changes which may subsequently increase the risk of disease transmission and weaken the disease control efforts.

The basic reproduction number of the model is

where $S_0$ and $M_0$ are the respective values of $S$ and $M$ at the disease-free equilibrium. The derivation of $\mathcal{R}_{0}$ is provided in Appendix A. We see that $\mathcal{R}_{0}$ includes two parts, representing the contributions from two respective transmission routes: one starts from the exposed, or asymptomatic infectious, individuals (denoted as $\mathcal{R}_{E0}$), and the other starts from the symptomatic infectious individuals (denoted as $\mathcal{R}_{I0}$). We refer to $\mathcal{R}_{E0}$ as the asymptomatic reproduction number, and $\mathcal{R}_{I0}$ as the symptomatic reproduction number.

We have also conducted an equilibrium analysis of system (2.1) by assuming that all the parameters are constants. The details are provided in Appendix B. In particular, the analysis indicates that in such a homogeneous setting, the disease will persist and become endemic in the long run if $\mathcal{R}_{0} > 1$. In practical applications, however, public health and economic policies typically change with time, and disease control strategies are adjusted accordingly. These variations, subsequently, impact the progression of the epidemic and lead to more complex situations than a homogeneous scenario. In the next section, we will use numerical simulation to explore such time-dependent behaviors of the system dynamics.

3.   Simulation results

We apply our model to the COVID-19 epidemic in the US state of Tennessee. We use the unemployment rate as the key economic indicator to infer the economic development level. Specifically, we set $C(t) = 1 - U(t)$ in our model, where $U(t)$ denotes the unemployment rate at time $t$. We have collected daily and weekly data for the COVID-19 cases and the unemployment rates in Tennessee ^[46,47], from March 28, 2020 to December 31, 2020. We fit our model to such data and conduct numerical simulation to examine the interaction between the epidemic spread, the disease management, and the economic development.

Since the epidemic progression of COVID-19 exhibits very different behaviors at different times, we divide the entire time interval (from March 28 to December 31, 2020) into five consecutive periods:

● Period 1: From March 28 to April 30

● Period 2: From May 1 to May 21

● Period 3: From May 22 to July 20

● Period 4: From July 21 to September 30

● Period 5: From October 1 to December 31

For each period, we conduct separate data fitting, with model parameters and initial conditions defined in Table 1. A gradient-based nonlinear constrained optimization procedure, implemented through the Matlab function fmincon, is employed to fit all the model parameters except the population influx rate $\lambda$, the recovery rate $\gamma$, the incubation period $\alpha^{-1}$, and the natural and disease-induced death rates $\mu$ and $\omega$, whose values are prescribed from the literature. We have found that this approach performs sufficiently well in our study. Nevertheless, we mention that many other data fitting techniques can be possibly applied, including some more sophisticated approaches such as a two-level fitting method presented in ^[52].

The fitting in different periods results in different estimates of parameter values in the model. All these fitted parameter values are listed in Table 2. We focus our attention on the data fitting and numerical simulation in Period 1, with details provided in Section 3.1. The analysis and discussion can be similarly applied to the other periods. A summary of the fitting and simulation results for Periods 2–5 is provided in Section 3.2.

3.1.   Fitting and simulation in Period 1

3.1.1.   Data fitting

Period 1, from March 28, 2020 to April 30, 2020, represents an initial phase of the COVID-19 epidemic in Tennessee (as well as in the entire US), with the number of infections fast increasing. Figure 2(a) shows our model fitting to the number of daily active cases reported by the Tennessee Department of Health ^[46], where we see that the number of active infections increased from below 1000 in the beginning to above 4000 at the end of the period. We observe a reasonably good match between our simulation result and the reported data.

Based on the estimated parameter values from Table 2, we evaluate the basic reproduction number from Eq (2.2):

where the subscript '1' indicates the association with Period 1. Clearly, $\mathcal{R}_{10} > 1$, which accounts for the fast increase of infections in this period. We note that the contribution from the symptomatic infectious individuals ($I$) is much higher than that from the asymptomatic infectious individuals ($E$) in shaping the disease risk in this period.

We plot the numerical curves of $S$, $E$, $I$, $R$, $M$ and $C$ in Figure 2(b)–(d). In particular, we observe that the value of the mitigation level $M$ quickly increases in the first five days (March 28 to April 2) and then remains steady afterwards. The initial increase of $M$ is obviously stimulated by the ascending phase of the epidemic. The stabilization afterwards is partly explained by the implementation of the state-wide Stay-at-Home order in Tennessee from April 1 to April 30. The Stay-at-Home order represents a high level of disease mitigation and outbreak management, and our model variable $M$ reaches and stays at this high level after April 2, as shown in Figure 2(d). Meanwhile, we note that even under the Stay-at-Home regulation, essential businesses (such as grocery stores, pharmacies, convenience stores, mail and shipping services, home repair, and automotive sales and repair) were still open in the US, and this mitigation policy is doubtlessly weaker than a complete lock-down. This may explain that the value of $M$ is stabilized at somewhere around $88\%$ but not higher.

In parallel with the fast initial increase of the mitigation level $M$, there is a decrease of the economic development level $C$ that starts from March 28 and continues until April 12, before it is finally stabilized. The economic development is negatively impacted by the disease control and management, and the reduction in the value of $C$ accompanies the growth of $M$. Figure 2(d) shows that the rate of decrease for $C$ (i.e., the downward slope) is smaller than the rate of increase for $M$ (i.e., the upward slope), indicating some sort of resistance of the economic system in response to adverse conditions. The extended reduction of $C$ beyond the increasing phase of $M$ represents a delayed effect of the high mitigation level on the economic development. The curve of $C$ is finally stabilized around $89.5\%$, consistent with the unemployment rate of $10.4\%$ reported in Tennessee for the last week of April ^[47].

3.1.2.   Special cases

A major difference between our model (2.1) and a conventional epidemic model is the incorporation of two variables, $M$ and $C$, related to economics. We now examine the impact of these two variables on the dynamics of system (2.1). To that end, we consider three special (and extreme) cases: $M = 0$, $C = 0$, and $M = C = 0$, with the same epidemic and economic datasets in Period 1 as described before.

For the first case, we remove the equation for $M$ from system (2.1), and set $M = 0$ in the remaining equations for $S$, $E$, $I$, $R$, and $C$. Similarly, for the second case, we remove the equation for $C$ from (2.1), and set $C = 0$ in the remaining equations for $S$, $E$, $I$, $R$, and $M$. For the third case, we remove the equations for $M$ and $C$ from (2.1), and set $M = C = 0$ in the remaining equations for $S$, $E$, $I$, and $R$. In each of these cases, we conduct data fitting for the corresponding reduced system and present the fitted parameter values in Table 3.

In order to compare the goodness-of-fit in different cases, we calculate the normalized root-mean-square error (NRMSE):

where $n$ is the number of days in this fitting period, $Y_{k}$ is the reported number of active cases on the $k$-th day, and $I_{k}$ is the numerical value of $I$ on the $k$-th day. The values of the NRMSE for the three special cases $M = 0$, $C = 0$, $M = C = 0$, and the original case $M, C\neq 0$, are found as 0.2036, 0.1287, 0.2218, and 0.1297, respectively. We see that, in particular, the first and third cases both lead to a lower quality of fitting than that associated with the original system (2.1).

We plot the result for the first case $M = 0$ in Figure 3. Comparing the parameter values in Table 3 to those in Table 2, we notice that the removal of the variable $M$ leads to a significant over-estimate of the parameter $d$, the natural reduction rate of the economic development. This results in a rapid (and unrealistic) decrease of $C$ to a level around $68\%$, as shown in Figure 3(b).

Meanwhile, we plot the result for the second case $C = 0$ in Figure 4. Comparing Table 3 to Table 2, we observe that the removal of the variable $C$ leads to an over-estimate of the parameter $p$, the natural reduction rate of the mitigation level, and an under-estimate of the parameter $b$, the implementation rate of the mitigation. These similarly result in a rapid (and unrealistic) decrease of $M$ to a level around $45\%$, as shown in Figure 4(b).

Figure 5 displays a histogram to compare the reproduction numbers in the three special cases ($M = 0$, $C = 0$, $M = C = 0$) and the original case with nonzero $M$ and $C$ for Period 1. In each case, we observe that the component $\mathcal{R}_{I0}$ is much higher than the other component $\mathcal{R}_{E0}$, showing that symptomatic infections contribute to most part of the disease risk in this period. Overall, the basic reproduction number $\mathcal{R}_{0}$ takes a lower value when $C = 0$ than that in the original case $M, C \neq 0$, indicating a lower disease risk in the extreme scenario where economic development is completely disregarded. On the other hand, $\mathcal{R}_{0}$ takes a higher value when $M = 0$ and $M = C = 0$ than that in the original case, indicating a higher disease risk in the extreme scenario where mitigation strategies are totally abandoned.

3.1.3.   Sensitivity analysis

To quantify the influence of model parameters on disease dynamics, we conduct a sensitivity analysis ^[53] for the parameters in system (2.1) using the data in Period 1. Table 4 displays the results for the relative sensitivity of the basic reproduction numbers, $\mathcal{R}_{10}$ in Period 1, in terms of model parameters. The relative sensitivity of $\mathcal{R}_{10}$ with respect to $\beta_E$, for example, is computed by

where Eq (2.2) has been used. Similar calculation applies to other parameters.

Figure 6 illustrates the variations of $\mathcal{R}_{10}\,$ in terms of the eight most sensitive parameters ranked in Table 4. Practically, among these eight parameters, at least five could be controlled through public health management: the transmission rates $\beta_I$ and $\beta_E$, the mitigation influx rate $\delta$, the mitigation decline rate $p$, and the mitigation implementation rate $b$. In particular, social distancing, mask wearing, and isolation of infected individuals would bring down the transmission rates $\beta_I$ and $\beta_E$ and could be most effective in reducing the disease risk. Meanwhile, strengthening epidemic management would increase the mitigation influx rate $\delta$ and implementation rate $b$, and promoting public awareness of the infection would weaken the mitigation decline rate $p$. These could also effectively reduce the value of the basic reproduction number. Additionally, improved treatment modalities, such as discovery of more effective therapeutic drugs for SARS-CoV-2, may be able to enlarge the recovery rate $\gamma$ (i.e., decrease the length of the average recovery period for COVID-19 patients) and help to reduce the overall risk of the epidemic.

Using the method described in ^[33], we have also computed the relative sensitivity of all the six state variables in model (2.1) with respect to these parameters, and the results are presented in Figure 7. In particular, the relative sensitivity of a state variable, $Y$, with respect to $\beta_E$ is given by $\frac{\partial Y}{\partial \beta_{E}}\cdot \frac{\beta_{E}}{Y} \,$. To obtain the partial derivative $\frac{\partial Y}{\partial \beta_{E}}$, we differentiate it with respect to $t$ to obtain

We then numerically solve for each $\frac{\partial Y}{\partial \beta_{E}}$ by connecting Eqs (2.1) and (3.4). The sensitivity computation for other parameters is similar.

We observe from Figure 7 that the variables $S$, $E$, $I$ and $R$ are all very sensitive to the parameters $\beta_I$, $b$, $\gamma$ and $\delta$. Meanwhile, all the state variables except $C$ have a high level of sensitivity with respect to $p$. Additionally, the variable $M$ is highly sensitive to $\delta$. On the other hand, $C$ is highly sensitive to $d$, the natural reduction rate of economic development, and $c_S$, the economic contribution rate from the susceptible individuals that would account for the major labor supply. These parameters, except $d$ and $c_S$, are also ranked among the most sensitive parameters for the basic reproduction number in Table 4. The fact that the variable $C$ and the basic reproduction number $\mathcal{R}_{10}$ have different sensitivity dependence on parameters is a reflection of the different focuses between economic development and public health mitigation.

3.1.4.   Simulated scenarios

Our sensitivity analysis in Section 3.1.3 has identified seven parameters: $\beta_I$, $b$, $\gamma$, $\delta$, $p$, $d$, and $c_S$, which are highly sensitive for the model outcome. We now conduct a detailed simulation to quantify the change of the model variables when the values of these parameters vary. We do not consider here the population influx rate $\lambda$ since it is a demographic characteristic that is largely independent of public health management and short-term economic development. We focus on the three state variables $I$, $M$ and $C$ in the discussion below. For comparison, we refer to the parameter values provided in Table 2 for Period 1 as their base values, and the simulation result presented in Figure 2 as the base scenario.

Variation of $\beta_{I}$. We first examine the impact of the symptomatic transmission rate $\beta_{I}$. We consider a scenario where $\beta_{I}$ is reduced to $10\%$ of its base value given in Table 2, and another scenario where $\beta_{I}$ is increased by $10$ times of its base value. Figure 8 displays the curves of $M$, $C$ and $I$ in these two hypothetical scenarios, where all the other parameters are fixed at values given in Table 2. The left panel of Figure 8 shows that, although the reduction of $\beta_{I}$ does not have much impact on $M$ and $C$, it substantially changes the behavior of $I$. As a result, the curve of $I$ keeps decreasing throughout Period 1 due to the significant reduction of disease transmission. On the other hand, the right panes shows that the 10-fold increase of $\beta_{I}$ leads to a dramatic increase in the number of infections, up to $3 \times 10^6$ which is more than $43\%$ of the total population in Tennessee. Consequently, this causes a sharp decrease in the economic development $C$, down to a level around $40 \%$, due to the significant loss of labor supply.

Variation of $\gamma$. Next, we consider the impact of the recovery rate $\gamma$. Figure 9 shows the simulation results when $\gamma$ is hypothetically reduced to $10\%$ (left panel), and hypothetically increased by $10$ times (right panel), in reference to its base value provided in Table 2. All other parameters are fixed. As can be naturally expected, the lower value of $\gamma$ leads to an increase of the number of active infections due to the prolonged disease recovery, while the higher value of $\gamma$ leads to a reduction of the number of active infections due to the shortened recovery period. In contrast, the variation of $\gamma$ has little influence on $M$ and $C$, indicating their low sensitivity on the recovery rate, which is consistent with the sensitivity results presented in Figure 7(e), (f).

Variation of $\delta$. We plot the simulation result with respect to the variation of the mitigation influx rate $\delta$ in Figure 10. As observed in Figure 7, $M$ and $I$ are positively and negatively correlated to $\delta$, respectively, and both $M$ and $I$ are highly sensitive to $\delta$, whereas $C$ has a low sensitivity level for $\delta$. The left panel of Figure 10 demonstrates that when $\delta$ is reduced to $10\%$ of its base value, $M$ rapidly decreases to a level around $9\%$, and $I$ quickly increases to a number about $5$ times of the base scenario in Figure 2c. The right panel illustrates that when $\delta$ is just increased by $10\%$ from its base value, $M$ would rise to a very high level (around $97\%$) and $I$ would be at a level lower than the base scenario in Figure 2c.

Variation of $c_S$. We learn from Figure 7(f) that $C$ is positively correlated and highly sensitive to the labor contribution rate $c_S$. Figure 11 illustrates the two hypothetical scenarios when $c_S$ is reduced to $10\%$ of its base value (left panel), where $C$ dramatically decreases to a level around $9\%$, and when $c_S$ is increased to $110\%$ of its base value (right panel), where $C$ reaches a very high level (around $98\%$). There is little change in the $M$ and $I$ curves when $c_S$ varies.

Variations of $p$, $d$ and $b$. Furthermore, we plot the simulation results for the variation of $p$ (the natural reduction rate of disease mitigation) in Figure 12, for the variation of $d$ (the natural reduction rate of economic development) in Figure 13, and for the variation of $b$ (the mitigation implementation rate) in Figure 14. The behaviors of $M$, $C$ and $I$ with respect to the variations of these parameters are consistent with the sensitivity results in Figure 7.

Additionally, we have combined the solution curves of $I$ in most of these simulated scenarios and presented the results in Figure 15, so that we may easily observe the quantitative impacts of these parameter variations on the number of active infections in Period 1. Among these presented scenarios, the reduction of $I$ is most significant with the (hypothetical) increase of the recovery rate $\gamma$ and the decrease of the transmission rate $\beta_I$. Note that a few off-the-scale scenarios are not included in this plot. Moreover, a histogram for the changes of the reproduction numbers under most of these parameter variations is presented in Figure 16, which provides a quantitative demonstration of the sensitivity prediction in Table 4. Consistent with the result in Figure 15, the decrease of $\beta_I$ and increase of $\gamma$ appear most effective in reducing the value of the basic reproduction number.

3.2.   Fitting and simulation in Periods 2–5

We have presented and discussed our detailed simulation results for Period 1 in Section 3.1. For completeness, we summarize the data fitting and numerical results for Periods 2–5 in Figures 17–20. These four periods represent the time from May 1, 2020 to December 31, 2020. The values of fitted parameters for each period are listed in Table 2. We observe several ups and downs for the curves of $I$, $M$ and $C$ throughout these periods, in accordance with the evolution of the infected cases and unemployment rates reported in Tennessee and the varied disease risk levels at different times. Based on the fitted parameters, the basic reproduction number for each of these periods is found as follows:

● Period 2: $\mathcal{R}_{20} = 0.22$

● Period 3: $\mathcal{R}_{30} = 1.69$

● Period 4: $\mathcal{R}_{40} = 0.91$

● Period 5: $\mathcal{R}_{50} = 1.38$

Comparing the basic reproduction numbers for all the five periods, we see that $\mathcal{R}_{10} = 2.01$ is the highest, which is associated with the initial ascending phase of the epidemic where the infections are typically increasing exponentially. The substantially reduced value of $\mathcal{R}_{20}$ for Period 2 is a result of the stay-at-home order in Tennessee that was implemented in April 2020, leading to a very low level of social contacts and disease transmission risk at the beginning of Period 2. As the stay-at-home order expired, businesses were re-opened, and people gradually resumed most of the their normal activities, our model shows that $\mathcal{R}_{30}$ bounces back to a high value for Period 3. This is followed by another reduction of the basic reproduction number for Period 4, where $\mathcal{R}_{40}$ is slightly below unity, possibly due to the improved public awareness and health management stimulated by the high disease risk in Period 3. Finally, the rebound of $\mathcal{R}_{50}$ for Period 5 is a reflection of the second COVID-19 wave that swept through the US during the last few months of 2020. Our simulation result for the evolution of the mitigation level $M$ shows a pattern consistent with that of the policy curves in Figure 1 (right).

In addition, we present the result of data fitting for the entire time interval (from March 28, 2020 to December 31, 2020) in Figure 21, where the values of fitted parameters are given by the last column in Table 2. We see that for several timeframes, some of which even span 1–2 months, our simulation curve is unable to well catch the behavior of the reported data, which provides a justification for our separated fitting and simulation in different periods.

4.   Discussion

We have presented a new mathematical model to study the interaction between the disease transmission and spread, the epidemic management, and the economic development associated with COVID-19. We have also fitted our model to the reported COVID-19 cases and unemployment rates in the US state of Tennessee to illustrate the application of our modeling framework. We have focused our attention on the model fitting and simulation during the initial period of the COVID-19 epidemic in Tennessee (from March 28 to April 30, 2020), though our simulation study for the remainder of the year 2020 (until December 31) is also summarized. We have considered a few special cases, including the extreme scenarios where there is no mitigation activity ($M = 0$) or there is no economic development ($C = 0$), and compared the different dynamics between those and the realistic scenario where both disease control and economic growth are present. Meanwhile, we have conducted a thorough sensitivity analysis of the model parameters, and utilized detailed simulation results to examine the impact on the model outcome when the values of several sensitive parameters are varied.

Our study represents an application of mathematical modeling and simulation at the interface between epidemiology and economics. The main contributions of this work can be summarized as follows.

Contribution to economics. Our findings could provide useful insight into economic development under the impact of COVID-19. Although our model is coarse-grained, simplifying the many factors in disease management and economic development to two variables ($M$ and $C$), the simulation results display a rich set of dynamics that highlight the complex interaction involved. Specifically, our numerical findings show several typical patterns: (i) increased disease prevalence stimulates stronger mitigation strategies; (ii) a higher level of disease control leads to reduced economic growth (measured by a higher unemployment rate); (iii) a lower level of economic development could help to slow down the epidemic spread, but the demand for economic recovery would tend to weaken the pandemic management and increase the risk of infection; (iv) within the triad of the epidemic progression, public health management, and economic development, the change of any single component will impact the other two, though such an impact may not be immediately noticeable and may often be seen in a delayed manner. These observations indicate the importance of striking a balance between the disease control and economic growth during the pandemic era, and could provide useful guidelines toward decision making and policy development.

Contribution to epidemic modeling. Our work represents a new modeling study in the epidemiology of COVID-19. Under a homogeneous setting, our mathematical analysis (presented in Appendix B) indicates that if $\mathcal{R}_0 > 1$, the disease would persist in the long term and the stable endemic equilibrium could potentially have a basin of attraction in the entire domain. Practically, however, public health management and change of human behavior could reduce the value of $\mathcal{R}_0$ below unity, so that the disease may be eradicated in a population. Even if this is not possible, strategical control measures could minimize the impact of the endemism, which, mathematically, corresponds to a decrease in the magnitude of the endemic equilibrium and/or a reduction of its basin of attraction to a small region. Furthermore, the interaction between economic growth and disease control may strongly impact the epidemic progression and lead to time-dependent system dynamics, so that the assumption of a homogeneous setting is no longer valid. Specifically, the number of reported COVID-19 cases in Tennessee (and possibly in the entire US as well) from March 2020 to December 2020 exhibits a few distinct features at different times, which has motivated us to divide the whole time interval (which is more than 9 months) into five different periods for data fitting and simulation. This allows us to better fit our model to the real data within each time frame. The simulation results at the end of each time period naturally provide initial conditions to the next period. Though we have focused on the first period, this separated fitting approach would enable us to conduct a detailed investigation of the epidemic-economic dynamics in each different period when needed, and the method can be easily extended to the fitting and simulation in other places and/or times.

Additionally, this modeling study builds a foundation for deeper predictive investigations into the relationship between the pandemic progression, disease control, and economic development. In particular, an interesting and practically meaningful task is to explore the 'best' balance between the pandemic management and economic growth. This can be formulated as an optimal control problem ^[54,55,56], and several epidemic-economic optimal policy studies have already been performed with various focuses ^{[57,58,59,60]}. For our modeling investigation, mathematical theory and numerical simulation can be applied to find a possible optimal control solution that could minimize the costs associated with the infection (morbidity and mortality), the control measures (prevention, tests, diagnosis, vaccination, treatment, etc.), and the reduction in economic development (increased unemployment rates, in particular). Such a result would provide helpful guidelines to the government, the public health administration, and other policy makers.

Acknowledgments

Research of JB is partially supported by Natural Science Youth Project of the Educational Department of Liaoning Province under grant number LQN202004, Research Project of Economic and Social Development in Liaoning Province under grant number 2022lslqnwzzkt-013, Young Women's Applied Mathematics Support Research Project of China Society for Industrial and Applied Mathematics. Research of XW is partially supported by Faculty Research and Creative Activity (RCA) Grant awarded by College of Arts and Sciences at the University of Tennessee at Chattanooga. Research of JW is partially supported by the National Institutes of Health under grant number 1R15GM131315.

Conflict of interest

The three authors, JB, XW and JW, declare that there is no conflict of interest in this work.

Appendix A.   Domain and reproduction number

We first determine a feasible domain for our mathematical model (2.1). If we add up the first four equations in system (2.1), we easily obtain

Meanwhile, adding the first two equations yields

Next, we derive conditions to ensure $0 \leq M \leq 1$ and $0 \leq C \leq 1$. For $M \geq 0$, we need $\frac{dM}{dt} \geq 0$ when $M = 0$. This yields $\delta +m I - f_{0}C \geq 0$, which holds for all $I \geq 0$ and $C \leq 1$ if

To ensure $M \leq 1$, we need $\frac{dM}{dt} \leq 0$ when $M = 1$. This leads to $\delta +m I -p -f_{0}C \leq 0$, which holds for all $I \leq \frac{\lambda}{\mu}$ and $C \geq 0$ if

With similar arguments, we can make certain $C \geq 0$ if

by using condition (A.2), and ensure $C \leq 1$ if

by using condition (A.1).

Thus, under conditions (A.3)–(A.6), the biologically feasible domain

is invariant with respect to the vector field of system (2.1).

Obviously, there is a unique disease-free equilibrium (DFE) of system (2.1) which is given by

Using conditions (A.3) and (A.6), we have $d\delta -f_{0}c_{S}S_{0} \geq f_0 (d - c_S \frac{\lambda}{\mu}) \geq f_0 (c_M \frac{\lambda}{\mu} - c_S \frac{\lambda}{\mu}) \geq 0$. Meanwhile, conditions (A.3) and (A.4) yield $p > \delta \geq f_0$, and conditions (A.5) and (A.6) yield $d \geq c_M \frac{\lambda}{\mu} > c_m \frac{\lambda}{\alpha +\mu} \geq g_0$. Hence, $dp -g_{0}f_{0} > 0$, and $pc_{S}S_{0}-g_{0} \delta > \delta (c_S S_0 - g_0) > \delta (c_m \frac{\lambda}{\alpha +\mu} - g_0) \geq 0$. These ensure that $M_0 \geq 0$ and $C_0 > 0$ in equation (A.8). Moreover, conditions (A.4) and (A.5) yield $d\delta -f_{0}c_{S}S_{0} < dp -g_{0}f_{0}$, or $M_0 < 1$, and conditions (A.3) and (A.6) yield $pc_{S}S_{0}-g_{0} \delta \leq dp -g_{0}f_{0}$, or $C_0 \leq 1$. Therefore, the DFE $x_0$ is within the domain $\Gamma$.

Using the DFE, we can evaluate the basic reproduction number of model (2.1) based on the standard next-generation matrix approach ^[61]. The new infection matrix $F$ and the transition matrix $V$ are

Thus, the next-generation matrix is given as

where $a_{11} = \frac{\beta_{E}S_{0}}{(1+bM_{0})(\alpha+\mu)}+\frac{\alpha\beta_{I}S_{0}}{(1+bM_{0}) (\alpha+\mu)(w+\gamma+\mu)}$ and $a_{12} = \frac{\beta_{I}S_{0}}{(1+bM_{0})(w+\gamma+\mu)}$.

The basic reproduction number is then computed by the spectral radius of the next-generation matrix: $\mathcal{R}_{0} = \rho(FV^{-1})$; i.e.,

Appendix B.   Equilibrium analysis

Here we analyze the equilibrium points of system (2.1) and their stability properties. We have determined the unique disease-free equilibrium in Appendix A. Now we establish the following result regarding the existence and uniqueness of the endemic equilibrium.

Theorem B.1. Assume $m \lambda - f_0 \gamma \geq 0$.When $\mathcal{R}_{0} > 1$, there exists a unique endemic equilibrium$x_{*} = (S_{*}, E_{*}, I_{*}, R_{*}, M_{*}, C_{*})$ defined by Eqs(B.7)–(B.12).

Proof. At an endemic equilibrium of system (2.1), we have

Adding (B.1) to (B.2), we obtain

By Eqs (B.3) and (B.4), we have

and

Using Eqs (B.5) and (B.6), we have

and

Substituting Eqs (B.7)–(B.11) to (B.2), we obtain a unique non-trivial solution for $E_*$:

where we may re-write $A$ as

We know $dp -f_{0}g_{0} > 0$ from Eq (A.8) and the analysis afterwards. Using conditions (A.3) and (A.6), we have $d \delta \frac{\alpha+\mu}{\lambda} - f_0 c_E \geq c_M \frac{\lambda}{\mu} \cdot \delta \frac{\alpha+\mu}{\lambda} - f_0 c_E > c_M \delta - f_0 c_E \geq 0$. Meanwhile, using the assumption $m \lambda - f_0 \gamma \geq 0$, we obtain $md - f_0 c_R \frac{\gamma}{\mu} \geq m c_M \frac{\lambda}{\mu} - f_0 c_R \frac{\gamma}{\mu} \geq \frac{c_R}{\mu} (m \lambda - f_0 \gamma) \geq 0$. Hence, $A > 0$, which yields $E_* > 0$ when $\mathcal{R}_{0} > 1$. Consequently, there is a unique endemic equilibrium $x_{*}$ when $\mathcal{R}_{0} > 1$.

Next, we consider the stability properties of the DFE and endemic equilibrium. Based on the standard theory of the basic reproduction number ^[61], the DFE is locally asymptotically stable when $\mathcal{R}_{0} < 1$ and becomes unstable when $\mathcal{R}_{0} > 1$. We proceed to study the stability of the endemic equilibrium near the bifurcation point $\mathcal{R}_{0} = 1$ using local bifurcation analysis ^[62]. We first introduce the following lemma.

Lemma B.1. (^[61], Theorem 4)Consider the disease transmission model defined by $\dot{x} = f(x, \mu)$, where$\mu$ is a bifurcation parameter such that $\mathcal{R}_{0} < 1$ for $\mu < 0$ and$\mathcal{R}_{0} > 1$ for $\mu > 0$ and such that $x_{0}$ is a DFE for all values of $\mu$.Assume that the zero eigenvalue of $D_{x}f(x_{0}, 0)$ is simple, where $D_{x}f(x_{0}, 0)$is the partial derivative of $f$ with respect to $x$ evaluated at the point $x = x_{0}$, $\mu = 0$, and let $v$ and $\omega$ be the corresponding left and right nullvectorschosen such that $v \omega = 1$.Let $a$ and $b$ be defined as follows:

and assume that $b \neq 0$. Also assume that the DFE is stable when $\mu < 0$.Then $b > 0$, and there exists $\varepsilon > 0$ such that

(i) if $a < 0$, there are locally asymptotically stable endemic equilibria near$x_{0}$ for $0 < \mu < \varepsilon$;

(ii) if $a > 0$, there are unstable endemic equilibria near $x_{0}$ for $-\varepsilon < \mu < 0$.

We utilize Lemma B.1 to prove the result below.

Theorem B.2. When $\mathcal{R}_{0} > 1$ and $\mathcal{R}_{0} -1$ is sufficiently small, the endemic equilibrium $x_{*}$ of system (2.1) is locally asymptotically stable.

Proof. Let us re-write system (2.1) as $\dot{x} = f(x, \mathcal{R}_{0})$, where $x = [ x_1, \, x_2, \, x_3, \, x_4, \, x_{5}, \, x_{6}]^T = [ S, \, E, \, I, \, R, \, M, \, C]^T$ and $f = [ f_1, \, f_2, \, f_3, \, f_4, \, f_{5}, \, f_{6} ]^T = [ \dot{S}, \, \dot{E}, \, \dot{I}, \, \dot{R}, \, \dot{M}, \, \dot{C} ]^T$.

The Jacobian matrix of system (2.1) at $(x_0, 1)$ is given by

Let $\mathcal{R}_{0} = 1$. We manipulate the characteristic polynomial as follows:

Note that $\mathcal{R}_{0}$ is defined in Eq (A.11). We define

Since $\mathcal{R}_{01} < \mathcal{R}_{0} = 1$, we have

Since $pd -f_{0}g_{0} > 0$, we see that the zero eigenvalue of $D_{x}f(x_{0}, 1)$ is simple and all other eigenvalues of $D_{x}f(x_{0}, 1)$ have negative real parts.

All the second derivatives of $f_{i}$ in (B.14) are zero at the DFE except the following: $\frac{\partial^{2} f_{1}}{\partial S \partial E}(x_0, 1) = \frac{-\beta_{E}}{1+bM_{0}}$, $\frac{\partial^{2} f_{1}}{\partial S \partial I}(x_0, 1) = \frac{-\beta_{I}}{1+bM_{0}}$, $\frac{\partial^{2} f_{1}}{\partial E \partial S}(x_0, 1) = \frac{-\beta_{E}}{1+bM_{0}}$, $\frac{\partial^{2} f_{1}}{\partial I \partial S}(x_0, 1) = \frac{-\beta_{I}}{1+bM_{0}}$, $\frac{\partial^{2} f_{1}}{\partial E \partial M}(x_0, 1) = \frac{\beta_{E}S_{0}b}{(1+bM_{0})^{2}}$, $\frac{\partial^{2} f_{1}}{\partial I \partial M}(x_0, 1) = \frac{\beta_{I}S_{0}b}{(1+bM_{0})^{2}}$, $\frac{\partial^{2} f_{1}}{\partial M \partial E}(x_0, 1) = \frac{\beta_{E}S_{0}b}{(1+bM_{0})^{2}}$, $\frac{\partial^{2} f_{1}}{\partial M \partial I}(x_0, 1) = \frac{\beta_{I}S_{0}b}{(1+bM_{0})^{2}}$, $\frac{\partial^{2} f_{2}}{\partial S \partial E}(x_0, 1) = \frac{\beta_{E}}{1+bM_{0}}$, $\frac{\partial^{2} f_{2}}{\partial S \partial I}(x_0, 1) = \frac{\beta_{I}}{1+bM_{0}}$, $\frac{\partial^{2} f_{2}}{\partial E \partial S}(x_0, 1) = \frac{\beta_{E}}{1+bM_{0}}$, $\frac{\partial^{2} f_{2}}{\partial I \partial S}(x_0, 1) = \frac{\beta_{I}}{1+bM_{0}}$, $\frac{\partial^{2} f_{2}}{\partial E \partial M}(x_0, 1) = \frac{-\beta_{E}S_{0}b}{(1+bM_{0})^{2}}$, $\frac{\partial^{2} f_{2}}{\partial I \partial M}(x_0, 1) = \frac{-\beta_{I}S_{0}b}{(1+bM_{0})^{2}}$, $\frac{\partial^{2} f_{2}}{\partial M \partial E}(x_0, 1) = \frac{-\beta_{E}S_{0}b}{(1+bM_{0})^{2}}$, $\frac{\partial^{2} f_{2}}{\partial M \partial I}(x_0, 1) = \frac{-\beta_{I}S_{0}b}{(1+bM_{0})^{2}}$.

Corresponding first derivatives of $f_{i}$ in (B.15) are listed as follows: $\frac{\partial f_1}{\partial E}(x_0, 1) = \frac{-\beta_{E}S_{0}}{1+bM_{0}} = -(\alpha+\mu)\mathcal{R}_{01} = -(\alpha+\mu)(\mathcal{R}_{0}-\mathcal{R}_{02})$, $\frac{\partial f_1}{\partial I}(x_0, 1) = \frac{-\beta_{I}S_{0}}{1+bM_{0}} = -\frac{(\alpha+\mu)(w+\gamma+\mu)}{\alpha}\mathcal{R}_{02} = -\frac{(\alpha+\mu)(w+\gamma+\mu)}{\alpha}(\mathcal{R}_{0}-\mathcal{R}_{01})$, $\frac{\partial f_2}{\partial E}(x_0, 1) = \frac{\beta_{E}S_{0}}{1+bM_{0}}-(\alpha+\mu) = (\alpha+\mu)(\mathcal{R}_{01}-1) = (\alpha+\mu)(\mathcal{R}_{0}-\mathcal{R}_{02}-1)$, $\frac{\partial f_2}{\partial I}(x_0, 1) = \frac{\beta_{I}S_{0}}{1+bM_{0}} = \frac{(\alpha+\mu)(w+\gamma+\mu)}{\alpha}\mathcal{R}_{02} = \frac{(\alpha+\mu)(w+\gamma+\mu)}{\alpha}(\mathcal{R}_{0}-\mathcal{R}_{01})$.

The non-zero second derivatives of $f_{i}$ in (B.15) are $\frac{\partial^{2} f_1}{\partial E\partial\mathcal{R}_{0}}(x_0, 1) = -(\alpha+\mu)$, $\frac{\partial^{2} f_1}{\partial I\partial\mathcal{R}_{0}}(x_0, 1) = -\frac{(\alpha+\mu)(w+\gamma+\mu)}{\alpha}$, $\frac{\partial^{2} f_2}{\partial E\partial\mathcal{R}_{0}}(x_0, 1) = \alpha+\mu$, $\frac{\partial^{2} f_2}{\partial I\partial\mathcal{R}_{0}}(x_0, 1) = \frac{(\alpha+\mu)(w+\gamma+\mu)}{\alpha}$.

We choose $v$ and $\omega$ such that they are orthogonal to $D_{x}f(x_{0}, 1)$ (i.e., $v\cdot D_{x}f(x_{0}, 1) = 0$, $D_{x}f(x_{0}, 1)\cdot\omega = 0$), and $v\cdot\omega = 1$. With some algebraic manipulation, we obtain $v = (0, v_{2}, v_{3}, 0, 0, 0)$, where $v_{2} = \frac{1}{1+\frac{\alpha\beta_{I}S_{0}}{1+bM_{0}}} > 0$, $v_{3} = \frac{\frac{\beta_{E}S_{0}}{1+bM_{0}}-(\alpha+\mu)}{-\alpha}v_{2} = \frac{\beta_{I}S_{0}}{(1+bM_{0})(w+\gamma+\mu)}v_{2}$ due to $\mathcal{R}_{0} = 1$; and $\omega = (\omega_{1}, \, \omega_{2}, \, \omega_{3}, \, \omega_{4}, \, \omega_{5}, \, \omega_{6})$, where $\omega_{1} = -\frac{\alpha+\mu}{\mu} < 0$, $\omega_{2} = 1$, $\omega_{3} = \frac{\alpha}{w+\gamma+\mu} > 0$, $\omega_{4} = \frac{\gamma\alpha}{\mu(w+\gamma+\mu)}$, $\omega_{5} = \frac{f_{0}c_{S}\frac{\alpha +\mu}{u}}{pd-f_{0}g_{0}}+\frac{md\cdot\frac{\alpha}{w+\gamma+\mu}-f_{0}c_{E}-f_{0}c_{R}\cdot \frac{\gamma}{\mu}\cdot\frac{\alpha}{w+\gamma+\mu}}{pd-f_{0}g_{0}} > 0$, $\omega_{6} = \frac{-pc_{S}\frac{\alpha+\mu}{\mu}-g_{0} \frac{m\alpha}{w+\gamma+\mu}+pc_{E}+pc_{R}\frac{\gamma\alpha} {\mu(w+\gamma+\mu)}}{pd-f_{0}g_{0}}$.

Hence, Eq (B.14) yields

Based on Eq (B.15), we can also verify that $b = v_{2}(\alpha+\mu)(\omega_{2}+\omega_{3}\frac{w+\gamma+\mu}{\alpha}) > 0$. Therefore, we conclude that when $\mathcal{R}_0-1$ changes from negative to positive, a positive and locally asymptotically stable equilibrium $x_*\,$ occurs.

Due to the high dimension and strong nonlinearity of system (2.1), we have not fully resolved the global stability of the DFE and the endemic equilibrium. Nevertheless, for a special case where the variable $M$ is fixed at any value between $0$ and $1$, we are able to prove the global asymptotic stability of the simplified system. For illustration, let us substitute $M = M_*$ into system (2.1). Consequently, the basic reproduction number for the reduced system is given by

where

Theorem B.3. Let $M = M_*$.When $\mathcal{R}_{0} < 1$, the DFE $x_0$ of the reduced systemfrom (2.1) is globally asymptotically stable.

Proof. When $M = M_{*}$, the second and third equations of (2.1) imply that

Let $Y = (E, \, I)$. Then system (B.19) yields

where the matrices $F_*$ and $V_*$ are given in Eq (B.18).

By the Perron-Frobenius theorem, there exists a positive left eigenvector $u$ of the positive matrix $V_*^{-1}F_*$ with respect to the eigenvalue $\mathcal{R}_{0} = \rho(F_*V_*^{-1}) = \rho(V_*^{-1}F_*)$. Define the Lyapunov function

Note that $L \geq 0$, and $L = 0$ if and only if $Y = 0$. Differentiating $L$ along the solution of the reduced system from (2.1), we obtain

If $\mathcal{R}_{0} < 1$, then $L' \leq 0$. Furthermore, $L' = 0$ leads to $u^T Y = 0$, which yields $E = I = 0$. Hence, the largest invariant set where $L' = 0$ is the DFE $x_0$. By LaSalle's Invariance Principle, the DFE is globally asymptotically stable.

Theorem B.4. Let $M = M_*$.When $\mathcal{R}_{0} > 1$, the unique endemic equilibrium $x_{*}$ of the reduced systemfrom (2.1) is globally asymptotically stable.

Proof. When $M = M_{*}$, we consider the first three equations of (2.1):

We introduce the following Lyapunov function:

It is easy to verify that $L_1 \geq 0$ and $L_1 = 0$ if and only if $(S, \, E, \, I) = (S_*, \, E_*, \, I_*)$. Meanwhile, the following equations hold:

Using Eqs (B.23)–(B.25), we calculate the derivative of $L_{1}$ along the solution of the subsystem (B.21) in what follows:

where those terms marked with the same type of symbols are canceled out. It is clear that $L_{1}'|_{(B.22)} = 0$ if and only if $(S, \, E, \, I) = (S_*, \, E_*, \, I_*)$. Hence, according to LaSalle's Invariance Principle, the endemic equilibrium $(S_{*}, \, E_{*}, \, I_{*})$ of subsystem (B.21) is globally asymptotically stable.

Through substitution of $(S_{*}, \, E_{*}, \, I_{*})$ into the equations for $R$ and $C$:

it is straightforward to observe that $(R_{*}, \, C_{*})$ of subsystem (B.26) is globally asymptotically stable. The proof is then completed by combining the results for subsystems (B.21) and (B.26).