Unreported Cases for Age Dependent COVID-19 Outbreak in Japan

Abstract

We investigate the age structured data for the COVID-19 outbreak in Japan. We consider a mathematical model for the epidemic with unreported infectious patient with and without age structure. In particular, we build a new mathematical model and a new computational method to fit the data by using age classes dependent exponential growth at the early stage of the epidemic. This allows to take into account differences in the response of patients to the disease according to their age. This model also allows for a heterogeneous response of the population to the social distancing measures taken by the local government. We fit this model to the observed data and obtain a snapshot of the effective transmissions occurring inside the population at different times, which indicates where and among whom the disease propagates after the start of public mitigation measures.

1. Introduction

COVID-19 disease caused by the severe acute respiratory syndrome coronavirus (SARS-CoV-2) first appeared in Wuhan, China, and the first cases were notified to WHO on 31 December 2019 [1,2]. Beginning in Wuhan as an epidemic, it then spread very quickly and was characterized a pandemic on 11 March 2020 [1]. Symptoms of this disease include fever, shortness of breath, cough, and a non-negligible proportion of infected individuals may develop severe forms of the symptoms leading to their transfer to intensive care units and, in some cases, death, see e.g., Guan et al. [3] and Wei et al. [4]. Both symptomatic and asymptomatic individuals can be infectious [4,5,6], which makes the control of the disease particularly challenging.

The virus is characterized by its rapid progression among individuals, most often exponential in the first phase, but also a marked heterogeneity in populations and geographic areas [7,8,9]. The number of reported cases worldwide exceeded 3 millions as of 3 May 2020 [10]. The heterogeneity of the number of cases and the severity according to the age groups, especially for children and elderly people, aroused the interest of several researchers [11,12,13,14,15]. Indeed, several studies have shown that the severity of the disease increases with the age and co-morbidity of hospitalized patients (see e.g., To et al. [15] and Zhou et al. [8]). Wu et al. [16] have shown that the risk of developing symptoms increases by $4\%$ per year in adults aged between 30 and 60 years old while Davies et al. [17] found that there is a strong correlation between chronological age and the likelihood of developing symptoms. Since completely asymptomatic individuals can also be contagious, a higher probability of developing symptoms does not necessarily imply greater infectiousness: Zou et al. [6] found that, in some cases, the viral load in asymptomatic patients was similar to that in symptomatic patients. Moreover while adults are more likely to develop symptoms, Jones et al. [18] found that the viral loads in infected children do not differ significantly from those of adults.

These findings suggest that a study of the dynamics of inter-generational spread is fundamental to better understand the spread of the coronavirus and most importantly to efficiently fight the COVID-19 pandemic. To this end the distribution of contacts between age groups in society (work, school, home, and other locations) is an important factor to take into account when modeling the spread of the epidemic. To account for these facts, some mathematical models have been developed [13,14,17,19,20]. In Ayoub et al. [19] the authors studied the dependence of the COVID-19 epidemic on the demographic structures in several countries but did not focus on the contacts distribution of the populations. In [13,14,17,20] a focus on the social contact patterns with respect to the chronological age has been made by using the contact matrices provided in Prem et al. [21]. While Ayoub et al. [19], Chikina and Pegden [20] and Davies et al. [17] included the example of Japan in their study, their approach is significantly different from ours. Indeed, Ayoub et al. [19] use a complex mathematical model to discuss the influence of the age structure on the infection in a variety of countries, mostly through the basic reproduction number ${\mathcal{R}}_0$. They use parameter values from the literature and from another study of the same group of authors [22], where the parameter identification is done by a nonlinear least-square minimization. Chikina and Pegden [20] use an age-structured model to investigate age-targeted mitigation strategies. They rely on parameter values from the literature and do discuss using age-structured temporal series to fit their model. Finally, Davies et al. [17] also discuss age-related effects in the control of the COVID epidemic, and use statistical inference to fit an age-structured SIR variant to data; the model is then used to discuss the efficiency of different control strategies. We provide a new, explicit computational solution for the parameter identification of an age-structured model. The model is based on the SIUR model developed in Liu et al. [23], which accounts for a differentiated infectiousness for reported and unreported cases (contrary to, for instance, other SIR-type models). In particular, our method is significantly different from nonlinear least-squares minimization and does not involve statistical inference.

In this article we focus on an epidemic model with unreported infectious symptomatic patients (i.e., with mild or no symptoms). Our goal is to investigate the age structured data of the COVID-19 outbreak in Japan. In Section 2 we present the age structured data and in Section 3 the mathematical models (with and without age structure). One of the difficulties in fitting the model to the data is that the growth rate of the epidemic is different in each age class, which lead us to adapt our early method presented in Liu et al. [23]. The new method is presented in the Appendix A. In Section 4 we present the comparison of the model with the data. In the last section we discuss our results.

2. Data

Patient data in Japan have been made public since the early stages of the epidemic with the quarantine of the Diamond Princess in the Haven of Yokohama. We used data from the website covid19japan.com (https://covid19japan.com. Accessed 6 May 2020) which is based on reports from national and regional authorities. Patients are labeled “confirmed” when tested positive to COVID-19 by PCR. Interestingly, the age class of the patient is provided for 13,660 out of 13,970 confirmed patients (97.8% of the confirmed population) as of 29 April. The age distribution of the infected population is represented in Figure 1 compared to the total population per age class (data from the Statistics Bureau of Japan estimate for 1 October 2019). In Figure 2 we plot the number of reported cases per 10,000 people of the same age class (i.e., the number of infected patients divided by the population of the age class times 10,000). Both datasets are given in

Table 1. The age distribution of Japan is taken from the Statistics Bureau of Japan [24]. The number of cases and the number of death the data come from Prefectural Governments and Japan Ministry of Health, Labour and Welfare.

Age group	$[0,10[$	$[10,20[$	$[20,30[$	$[30,40[$	$[40,50[$	$[50,60[$	$[60,70[$	$[70,80[$	$[80,90[$	$[90,100[$
Age class for 2019	9,859,515	11,171,044	12,627,964	14,303,042	18,519,755	16,277,853	16,231,582	15,926,926	8,939,954	2,309,313
Age class per 10,000 people	781	885	1000	1133	1467	1290	1286	1262	709	183
Confirmed Cases	211	327	2216	2034	2220	2355	1566	1289	857	304
Death	0	0	0	2	6	4	7	37	49	9

and a statistical summary is provided by

Table 2. Statistical summary of the data from Table 1.

Dataset	Japanese Population	Infected	Deceased
First Quartile	28	28	68
Median	48	44	75
Third Quartile	67	59	81

. Note that the high proportion of 20–60 years old confirmed patients may indicate that the severity of the disease is lower for those age classes than for older patients, and therefore the disease transmits more easily in those age classes because of a higher number of asymptomatic individuals. Elderly infected individuals might transmit less because they are identified more easily. The cumulative number of death (Figure 3) is another argument in favor of this explanation. We also reconstructed the time evolution of the reported cases in Figure 4 and Figure 5. Note that the steepest curves precisely concern the 20–60-year old, probably because they are economically active and therefore have a high contact rate with the population.

3. Methods

3.1. SIUR Model

The model consists of the following system of ordinary differential equations:

$$\left\{\begin{array}{c}{S}^{\prime }\left(t\right)=-\tau \left(t\right)S\left(t\right){\displaystyle \frac{{\displaystyle I\left(t\right)+U\left(t\right)}}{{\displaystyle N}}},\hfill \\ I^{\prime }\left(t\right)=\tau \left(t\right)S\left(t\right){\displaystyle \frac{{\displaystyle I\left(t\right)+U\left(t\right)}}{{\displaystyle N}}}-\nu I\left(t\right),\hfill \\ R^{\prime }\left(t\right)={\nu }_{1}I\left(t\right)-\eta R\left(t\right),\hfill \\ U^{\prime }\left(t\right)={\nu }_{2}I\left(t\right)-\eta U\left(t\right).\hfill \end{array}\right.$$

(1)

This system is supplemented by initial data

$$S\left(t_0\right)=S_0\ge 0,I\left(t_0\right)=I_0\ge 0,R\left(t_0\right)\ge 0 \mathrm{and} U\left(t_0\right)=U_0\ge 0.$$

(2)

Here $t\ge t_0$ is time in days, $t_0$ is the starting date of the epidemic in the model, $S\left(t\right)$ is the number of individuals susceptible to infection at time t, $I\left(t\right)$ is the number of asymptomatic infectious individuals at time t, $R\left(t\right)$ is the number of reported symptomatic infectious individuals at time t, and $U\left(t\right)$ is the number of unreported symptomatic infectious individuals at time t. A flow chart of the model is presented in Figure 6.

Asymptomatic infectious individuals $I\left(t\right)$ are infectious for an average period of $1/\nu$ days. Reported symptomatic individuals $R\left(t\right)$ are infectious for an average period of $1/\eta$ days, as are unreported symptomatic individuals $U\left(t\right)$. We assume that reported symptomatic infectious individuals $R\left(t\right)$ are reported and isolated immediately, and cause no further infections. The asymptomatic individuals $I\left(t\right)$ can also be viewed as having a low-level symptomatic state. All infections are acquired from either $I\left(t\right)$ or $U\left(t\right)$ individuals. A summary of the parameters involved in the model is presented in

Table 3. Parameters of the model.

Symbol	Interpretation	Method
$t_0$	Time at which the epidemic started	fitted
$S_0$	Number of susceptible at time $t_0$	fixed
$I_0$	Number of asymptomatic infectious at time $t_0$	fitted
$U_0$	Number of unreported symptomatic infectious at time $t_0$	fitted
$\tau \left(t\right)$	Transmission rate at time t	fitted
D	First day of public intervention	fitted
$\mu$	Intensity of the public intervention	fitted
$1/\nu$	Average time during which asymptomatic infectious are asymptomatic	fixed
f	Fraction of asymptomatic infectious that become reported symptomatic infectious	fixed
${\nu }_1=f\nu$	Rate at which asymptomatic infectious become reported symptomatic	fixed
${\nu }_2=(1-f)\nu$	Rate at which asymptomatic infectious become unreported symptomatic	fixed
$1/\eta$	Average time symptomatic infectious have symptoms	fixed

.

Our study begins in the second phase of the epidemics, i.e., after the pathogen has succeeded in surviving in the population. During this second phase $\tau \left(t\right)\equiv {\tau }_0$ is constant. When strong government measures such as isolation, quarantine, and public closings are implemented, the third phase begins. The actual effects of these measures are complex, and we use a time-dependent decreasing transmission rate $\tau \left(t\right)$ to incorporate these effects. The formula for $\tau \left(t\right)$ is

$$\left\{\begin{array}{c}\tau \left(t\right)={\tau }_0,0\le t\le D,\hfill \\ \tau \left(t\right)={\tau }_{0}exp\left(-\mu \left(t-D\right)\right),D<t.\hfill \end{array}\right.$$

(3)

The date D is the first day of public intervention and $\mu$ characterises the intensity of the public intervention.

A similar model has been used to describe the epidemics in mainland China, South Korea, Italy, and other countries, and give reasonable trajectories for the evolution of the epidemic based on actual data [23,25,26,27,28,29]. Compared with these models, we added a scaling with respect to the total population size N, for consistency with the age-structured model (12). This only changes the value of the parameter $\tau$ and does not impact the qualitative or quantitative behavior of the model.

3.2. Comparison of the Model (1) with the Data

At the early stages of the epidemic, the infectious components of the model $I\left(t\right)$, $U\left(t\right)$ and $R\left(t\right)$ must be exponentially growing. Therefore, we can assume that

$$I\left(t\right)=I_{0}exp\left({\chi }_2\left(t-t_0\right)\right).$$

The cumulative number of reported symptomatic infectious cases at time t, denoted by $CR\left(t\right)$, is

$${\displaystyle CR\left(t\right)={\nu }_1\underset{t_0}{\overset{t}{\int }}I\left(s\right)ds.}$$

(4)

Since $I\left(t\right)$ is an exponential function and $CR\left(t_0\right)=0$ it is natural to assume that $CR\left(t\right)$ has the following special form:

$$CR\left(t\right)={\chi }_{1}exp\left({\chi }_{2}t\right)-{\chi }_3.$$

(5)

As in our early articles [23,26,27,28,29], we fix ${\chi }_3=1$ and we evaluate the parameters ${\chi }_1$ and ${\chi }_2$ by using an exponential fit to

$${\chi }_{1}exp\left({\chi }_{2}t\right)\simeq C{R}_{data}\left(t\right).$$

We use only early data for this part, from day $t=d_1$ until day $t=d_2$, because we want to catch the exponential growth of the early epidemic and avoid the influence of saturation arising at later stages.

Remark 1.

The estimated parameters ${\chi }_1$ and ${\chi }_2$ will vary if we change the interval $\left[d_1,d_2\right]$.

Once ${\chi }_1,{\chi }_2,{\chi }_3$ are known, we can compute the starting time of the epidemic $t_0$ from (5) as:

$$CR\left(t_0\right)=0\iff {\chi }_{1}exp\left({\chi }_2{t}_0\right)-{\chi }_3=0\Rightarrow t_0={\displaystyle \frac{1}{{\chi }_2}}\left(ln\left({\chi }_3\right)-ln\left({\chi }_1\right)\right).$$

We fix $S_0=126.8\times {10}^6$, which corresponds to the total population of Japan. The quantities $I_0$, $R_0$, and $U_0$ correspond to the values taken by $I\left(t\right)$, $R\left(t\right)$ and $U\left(t\right)$ at $t=t_0$ (and in particular $R_0$ should not be confused with the basic reproduction number ${\mathcal{R}}_0$). We fix the fraction f of symptomatic infectious cases that are reported. We assume that between $80\%$ and $100\%$ of infectious cases are reported. Thus, f varies between $0.8$ and 1. We assume that the average time during which the patients are asymptomatic infectious $1/\nu$ varies between 1 day and 7 days. We assume that the average time during which a patient is symptomatic infectious $1/\eta$ varies between 1 day and 7 days. In other words we fix the parameters f, $\nu$, $\eta$. Since f and $\nu$ are known, we can compute

$${\nu }_1=f\nu \mathrm{and}{\nu }_2=(1-f)\nu .$$

(6)

Computing further (see below for more details), we should have

$$I_0={\displaystyle \frac{{\chi }_1{\chi }_{2}exp\left({\chi }_2{t}_0\right)}{f\nu }}={\displaystyle \frac{{\chi }_3{\chi }_2}{f\nu }},$$

(7)

$$\tau =N\frac{{\chi }_2+\nu }{S_0}{\displaystyle \frac{\eta +{\chi }_2}{{\nu }_2+\eta +{\chi }_2}},$$

(8)

$$R_0={\displaystyle \frac{{\nu }_1}{\eta +{\chi }_2}}{I}_0={\displaystyle \frac{f\nu }{\eta +{\chi }_2}}{I}_0,$$

(9)

and

$$U_0={\displaystyle \frac{{\nu }_2}{\eta +{\chi }_2}}{I}_0={\displaystyle \frac{(1-f)\nu }{\eta +{\chi }_2}}{I}_0.$$

(10)

By using the approach described in Diekmann et al. [30], van den Driessche and Watmough [31], the basic reproductive number for model (1) is given by

$${\mathcal{R}}_0={\displaystyle \frac{\tau S_0}{\nu N}}\left(1+{\displaystyle \frac{{\nu }_2}{\eta }}\right).$$

By using (8) we obtain

$${\mathcal{R}}_0=\frac{{\chi }_2+\nu }{\nu }{\displaystyle \frac{(\eta +{\chi }_2)}{{\nu }_2+\eta +{\chi }_2}}\left(1+{\displaystyle \frac{{\nu }_2}{\eta }}\right).$$

(11)

3.3. Model SIUR with Age Structure

In what follows we will denote $N_1,\dots ,N_{10}$ the number of individuals respectively for the age classes $[0,10[,\dots ,[90,100[$. The model for the number of susceptible individuals $S_1\left(t\right),\dots ,S_{10}\left(t\right)$, respectively for the age classes $[0,10[,\dots ,[90,100[$, is the following

$$\left\{\begin{array}{c}{S}_1^{\prime }\left(t\right)=-{\tau }_1{S}_1\left(t\right)\left[{\varphi }_{1,1}{\displaystyle \frac{\left(I_1\left(t\right)+U_1\left(t\right)\right)}{N_1}}+\dots +{\varphi }_{1,10}{\displaystyle \frac{\left(I_{10}\left(t\right)+U_{10}\left(t\right)\right)}{N_{10}}}\right],\hfill \\ ⋮\hfill \\ S_{10}^{\prime }\left(t\right)=-{\tau }_{10}{S}_{10}\left(t\right)\left[{\varphi }_{10,1}{\displaystyle \frac{\left(I_1\left(t\right)+U_1\left(t\right)\right)}{N_1}}+\dots +{\varphi }_{10,10}{\displaystyle \frac{\left(I_{10}\left(t\right)+U_{10}\left(t\right)\right)}{N_{10}}}\right].\hfill \end{array}\right.$$

(12)

The model for the number of asymptomatic infectious individuals $I_1\left(t\right),\dots ,I_{10}\left(t\right)$, respectively for the age classes $[0,10[,\dots ,[90,100[$, is the following

$$\left\{\begin{array}{c}{I}_1^{\prime }\left(t\right)={\tau }_1{S}_1\left(t\right)\left[{\varphi }_{1,1}{\displaystyle \frac{\left(I_1\left(t\right)+U_1\left(t\right)\right)}{N_1}}+\dots +{\varphi }_{1,10}{\displaystyle \frac{\left(I_{10}\left(t\right)+U_{10}\left(t\right)\right)}{N_{10}}}\right]-\nu I_1\left(t\right),\hfill \\ ⋮\hfill \\ I_{10}^{\prime }\left(t\right)={\tau }_{10}{S}_{10}\left(t\right)\left[{\varphi }_{10,1}{\displaystyle \frac{\left(I_1\left(t\right)+U_1\left(t\right)\right)}{N_1}}+\dots +{\varphi }_{10,10}{\displaystyle \frac{\left(I_{10}\left(t\right)+U_{10}\left(t\right)\right)}{N_{10}}}\right]-\nu I_{10}\left(t\right).\hfill \end{array}\right.$$

(13)

The model for the number of reported symptomatic infectious individuals $R_1\left(t\right),\dots ,R_{10}\left(t\right)$, respectively for the age classes $[0,10[,\dots ,[90,100[$, is

$$\left\{\begin{array}{c}{R}_1^{\prime }\left(t\right)={\nu }_1^1{I}_1\left(t\right)-\eta R_1\left(t\right),\hfill \\ ⋮\hfill \\ R_{10}^{\prime }\left(t\right)={\nu }_1^{10}{I}_{10}\left(t\right)-\eta R_{10}\left(t\right).\hfill \end{array}\right.$$

(14)

Finally the model for the number of unreported symptomatic infectious individuals $U_1\left(t\right),\dots ,U_{10}\left(t\right)$, respectively in the age classes $[0,10[,\dots ,[90,100[$, is the following

$$\left\{\begin{array}{c}{U}_1^{\prime }\left(t\right)={\nu }_2^1{I}_1\left(t\right)-\eta U_1\left(t\right),\hfill \\ ⋮\hfill \\ U_{10}^{\prime }\left(t\right)={\nu }_2^{10}{I}_{10}\left(t\right)-\eta U_{10}\left(t\right).\hfill \end{array}\right.$$

(15)

In each age class $[0,10[,\dots ,[90,100[$ we assume that there is a fraction $f_1,\dots ,f_{10}$ of asymptomatic infectious individual who become reported symptomatic infectious (i.e., with severe symptoms) and a fraction $(1-f_1),\dots ,(1-f_{10})$ who become unreported symptomatic infectious (i.e., with mild symptoms). Therefore we define

$$\begin{array}{c}{\nu }_1^1=\nu f_1\mathrm{and}{\nu }_2^1=\nu (1-f_1),\\ ⋮\\ {\nu }_1^{10}=\nu f_{10}\mathrm{and}{\nu }_2^{10}=\nu (1-f_{10}).\end{array}$$

(16)

In this model ${\tau }_1,\dots ,{\tau }_{10}$ are the respective transmission rates for the age classes $[0,10[,\dots ,[90,100[$.

The matrix ${\varphi }_{ij}$ represents the probability for an individual in the class i to meet an individual in the class j. In their survey, Prem and co-authors [21] present a way to reconstruct contact matrices from existing data and provide such contact matrices for a number of countries including Japan. Based on the data provided by Prem et al. [21] for Japan we construct the contact probability matrix $\varphi$. More precisely, we inferred contact data for the missing age classes $[80,90[$ and $[90,100[$. The precise method used to construct the contact matrix $\gamma$ is detailed in Appendix B. An analogous contact matrix for Japan has been proposed by Munasinghe, Asai and Nishiura [32]. The contact matrix $\gamma$ we used is the following

$$\begin{array}{c}\left[{\gamma }_{ij}\right]=\left[\begin{array}{cccccccccc}4.03& 0.92& 0.47& 1.69& 0.83& 0.92& 0.78& 0.56& 0.57& 0.57\\ 0.71& 8.06& 1.38& 1.36& 1.96& 1.74& 0.75& 0.86& 0.74& 0.57\\ 0.55& 1.05& 4.63& 2.25& 1.84& 1.92& 0.94& 0.46& 0.74& 0.73\\ 1.52& 1.20& 2.54& 4.97& 2.98& 2.40& 1.76& 0.99& 0.53& 0.73\\ 0.69& 1.42& 1.93& 2.87& 3.91& 2.76& 1.35& 1.33& 0.95& 0.53\\ 0.34& 0.48& 1.20& 1.46& 1.61& 2.97& 1.40& 0.98& 1.23& 0.95\\ 0.28& 0.18& 0.20& 0.52& 0.38& 0.77& 2.67& 1.72& 0.92& 1.23\\ 0.12& 0.10& 0.09& 0.18& 0.19& 0.25& 0.76& 1.99& 1.18& 0.93\\ 0.09& 0.10& 0.08& 0.09& 0.13& 0.17& 0.27& 0.64& 1.61& 1.19\\ 0.09& 0.09& 0.10& 0.08& 0.09& 0.13& 0.17& 0.27& 0.64& 1.61\end{array}\right],\end{array}$$

(17)

where the ith line of the matrix ${\gamma }_{ij}$ is the average number of contact made by an individuals in the age class i with an individual in the age class j during one day. Notice that the higher number of contacts are achieved within the same age class. The matrix of conditional probability $\varphi$ of contact between age classes is given by (18) and we plot a visual representation of this matrix in Figure 7.

$$\begin{array}{c}\left[{\varphi }_{ij}\right]=\left[\begin{array}{cccccccccc}0.35& 0.08& 0.04& 0.14& 0.07& 0.08& 0.06& 0.04& 0.05& 0.05\\ 0.03& 0.44& 0.07& 0.07& 0.10& 0.09& 0.04& 0.04& 0.04& 0.03\\ 0.03& 0.06& 0.30& 0.14& 0.12& 0.12& 0.06& 0.03& 0.04& 0.04\\ 0.07& 0.06& 0.12& 0.25& 0.15& 0.12& 0.08& 0.05& 0.02& 0.03\\ 0.03& 0.07& 0.10& 0.16& 0.22& 0.15& 0.07& 0.07& 0.05& 0.03\\ 0.02& 0.03& 0.09& 0.11& 0.12& 0.23& 0.11& 0.07& 0.09& 0.07\\ 0.03& 0.02& 0.02& 0.05& 0.04& 0.08& 0.30& 0.19& 0.10& 0.13\\ 0.02& 0.01& 0.01& 0.03& 0.03& 0.04& 0.13& 0.34& 0.20& 0.16\\ 0.02& 0.02& 0.01& 0.02& 0.02& 0.03& 0.06& 0.14& 0.36& 0.27\\ 0.02& 0.02& 0.03& 0.02& 0.02& 0.03& 0.05& 0.08& 0.19& 0.48\end{array}\right].\end{array}$$

(18)

4. Results

4.1. Model without Age Structure

The daily number of reported cases from the model can be obtained by computing the solution of the following equation:

$$D{R}^{\prime }\left(t\right)={\nu }_{1}I\left(t\right)-DR\left(t\right),\mathrm{for}t\ge t_0\mathrm{and}DR\left(t_0\right)=D{R}_0.$$

(19)

In Figure 8 and Figure 9 we employ the method presented previously in Liu et al. [29] to fit the data for Japan without age structure.

The model to compute the cumulative number of death from the reported individuals is the following

$$D^{\prime }\left(t\right)={\eta }_{D}pR\left(t\right),\mathrm{for}t\ge t_0\mathrm{and}D\left(t_0\right)=0,$$

(20)

where ${\eta }_D$ is the death rate of reported infectious symptomatic individuals and p is the case fatality rate (namely the fraction of death per reported infectious individuals).

In the simulation we chose $1/{\eta }_D=6$ days and the case fatality rate $p=0.286$ is computed by using the cumulative number of confirmed cases and the cumulative number of deaths (as of 29 April) as follows

$$p=\frac{\mathrm{cumulative}\mathrm{number}\mathrm{of}\mathrm{deaths}}{\mathrm{cumulative}\mathrm{number}\mathrm{of}\mathrm{reported}\mathrm{cases}}=\frac{393}{13744}.$$

(21)

In Figure 10 we plot the cumulative number of $D\left(t\right)$ by using the same simulations than in Figure 8 and Figure 9.

4.2. Model with Age Structure

In order to describe the confinement for the age structured model (12)–(15) we will use for each age class $i=1,\dots ,10$ a different transmission rate having the following form

$$\left\{\begin{array}{c}{\tau }_i\left(t\right)={\tau }_i,0\le t\le D_i,\hfill \\ {\tau }_i\left(t\right)={\tau }_{i}exp\left(-{\mu }_i\left(t-D_i\right)\right),D_i<t.\hfill \end{array}\right.$$

(22)

The date $D_i$ is the first day of public intervention for the age class i and ${\mu }_i$ is the intensity of the public intervention for each age class.

In Figure 11 we plot the cumulative number of reported cases as given by our model (12)–(15) (solid lines), compared with reported cases data (black dots). We used the method described in the Appendix A to estimate the parameters ${\tau }_i$ from the data. In Figure 12 we plot the cumulative number of unreported cases (solid lines) as given by our model with the same parameter values, compared to the existing data of reported cases (black dots).

In order to understand the role of transmission network between age groups in this epidemic, we plot in Figure 13 the transmission matrices computed at different times. The transmission matrix is the following

$$C\left(t\right)=\mathrm{diag}\left({\tau }_1\left(t\right),{\tau }_2\left(t\right),\dots ,{\tau }_{10}\left(t\right)\right)\times \varphi ,$$

(23)

where the matrix $\varphi$ describes contacts and is given in (18), and the transmission rates are the ones fitted to the data as in Figure 11

$${\tau }_i\left(t\right)={\tau }_i^0\left(t\right)exp(-{\mu }_i{(t-D_i)}_{+}).$$

During the early stages of the epidemic, the transmission seems to be evenly distributed among age classes, with a little bias towards younger age classes (Figure 13a). Younger age classes seem to react more quickly to social distancing policies than older classes, therefore their transmission rate drops rapidly (Figure 13b,c); one month after the start of social distancing measures, the transmission mostly occurs within elderly classes (60–100 years, Figure 13d).

5. Discussion

The recent COVID-19 pandemic has lead many local governments to enforce drastic control measures in an effort to stop its progression. Those control measures were often taken in a state of emergency and without any real visibility concerning the later development of the epidemics, to prevent the collapse of the health systems under the pressure of severe cases. Mathematical models can precisely help see more clearly what could be the future of the pandemic provided that the particularities of the pathogen under consideration are correctly identified. In the case of COVID-19, one of the features of the pathogen which makes it particularly dangerous is the existence of a high contingent of unidentified infectious individuals who spread the disease without notice. This makes non-intensive containment strategies such as quarantine and contact-tracing relatively inefficient but also renders predictions by mathematical models particularly challenging.

Early attempts to reconstruct the epidemics by using SIUR models were performed in Liu et al. [23,26,27,28], who used them to fit the behavior of the epidemics in many countries, by including undetected cases into the mathematical model. Here we extend our modeling effort by adding the time series of deaths into the equation. In Section 4 we present an additional fit of the number of disease-induced deaths coming from symptomatic (reported) individuals (see Figure 10). In order to fit properly the data, we were forced to reduce the length of stay in the R-compartment to 6 days (on average), meaning that death induced by the disease should occur on average faster than recovery. A shorter period between infection and death (compared to remission) has also been observed, for instance, by Verity et al. [7].

The major improvement in this article is to combine our early SIUR model with chronological age. Early results using age structured SIR models were obtained by Kucharski et al. [33] but no unreported individuals were considered and no comparison with age-structured data were performed. Indeed in this article we provide a new method to fit the data and the model. The method extends our previous method for the SIUR model without age (see Appendix A).

The data presented in Section 2 suggests that the chronological age plays a very important role in the expression of the symptoms. The largest part of the reported patients are between 20 and 60 years old (see Figure 1), while the largest part of the deceased are between 60 and 90 years old (see Figure 3). This suggests that the symptoms associated with COVID-19 infection are more severe in elderly patients, which has been reported in the literature several times (see e.g., Lu et al. [12], Zhou et al. [8]). In particular, the probability of being asymptomatic (our parameter f) should in fact depend on the age class.

Indeed, the best match for our model (see Figure 11) was obtained under the assumption that the proportion of symptomatic individual among the infected increases with the age of the patient. This linear dependency of f as a function of age is consistent with the observations of Wu et al. [16] that the severity of the symptoms increase linearly with age. As a consequence, unreported cases are a majority for young age classes (for age classes less than 50 years) and become a minority for older age classes (more than 50 years), see Figure 12. Moreover, our model reveals the fact that the policies used by the government to reduce contacts between individuals have strongly heterogeneous effects depending on the age classes. Plotting the transmission matrix at different times (see Figure 13) shows that younger age classes react more quickly and more efficiently than older classes. This may be due to the fact that the number of contacts in a typical day is higher among younger individuals. As a consequence, we predict that one month after the effective start of public measures, the new transmissions will almost exclusively occur in elderly classes. The observation that younger ages classes play a major roles in the transmission of the disease has been highlighted several times in the literature, see e.g., Davies et al. [17], Cao et al. [11], Kucharski et al. [33] for the COVID-19 epidemic, but also Mossong et al. [34] in a more general context.

We develop a new model for age-structured epidemic and provided a new and efficient method to identify the parameters of this model based on observed data. Our method differs significantly from the existing nonlinear least-squares and statistical inference methods and we believe that it produces high-quality results. Moreover, we only use the initial phase of the epidemic for the identification of the epidemiological parameters, which shows that the model itself is consistent with the observed phenomenon and argues against overfitting. Yet our study could be improved in several direction. We only use reported cases which were confirmed by PCR tests, and therefore the number of tests performed could introduce a bias in the observed data – and therefore our results. We are currently working on an integration of this number of tests in our model. We use a phenomenological model to describe the response of the population in terms of number of contacts to the mitigation measures imposed by the government. This could probably be described more precisely by investigating the mitigation strategies in terms of social network. Nevertheless we believe that our study offers a precise and robust mathematical method which adds to the existing literature.