A Hybrid Modeling Technique of Epidemic Outbreaks with Application to COVID-19 Dynamics in West Africa

Simple Summary

The intrinsic dynamics of the propagation of a disease changes along an epidemic course, especially for long lasting epidemics such as the COVID-19. Indeed, the natural evolution of the pathogen and countermeasures such as quarantining, lockdown, social distancing and vaccination modify the transmission dynamics of the disease. With a view to match these theoretical changes to potential changes in observed epidemiological data, we designed a hybrid modeling framework where we integrated: (1) two growth curves for daily reported positive cases, differentiating the early epidemic phase and a second phase with a potentially different dynamics; (2) two logistic regression models for daily recoveries and deaths; and (3) a SIQR (Susceptible, Infective, Quarantined, Recovered) mechanistic model to provide an overview of the dynamics of the disease in the target population. This joint modeling approach allows explicit analytical expressions for the different compartments of the SIQR model, circumventing common identifiability issues in such models. The changes in the disease transmission pattern can be subjected to countermeasures so as to assess their effectiveness along time. For illustrative purposes, we applied the approach to COVID-19 data from West Africa. It turned out that the first imported COVID-19 case(s) in West Africa likely entered the region between 28 January and 7 February 2020. Moreover, the first measures implemented by West African authorities impacted the dynamics of the disease one month after the outbreak.

Abstract

The widely used logistic model for epidemic case reporting data may be either restrictive or unrealistic in presence of containment measures when implemented after an epidemic outbreak. For flexibility in epidemic case reporting data modeling, we combined an exponential growth curve for the early epidemic phase with a flexible growth curve to account for the potential change in growth pattern after implementation of containment measures. We also fitted logistic regression models to recoveries and deaths from the confirmed positive cases. In addition, the growth curves were integrated into a SIQR (Susceptible, Infective, Quarantined, Recovered) model framework to provide an overview on the modeled epidemic wave. We focused on the estimation of: (1) the delay between the appearance of the first infectious case in the population and the outbreak (“epidemic latency period”); (2) the duration of the exponential growth phase; (3) the basic and the time-varying reproduction numbers; and (4) the peaks (time and size) in confirmed positive cases, active cases and new infections. The application of this approach to COVID-19 data from West Africa allowed discussion on the effectiveness of some containment measures implemented across the region.

1. Introduction

The ravages of the COVID-19 pandemic has deepened the need for mathematical and statistical tools to understand the dynamics of epidemics across the world. Simple mathematical models of infectious diseases are useful for providing insight into epidemic trajectories and disease dynamics [1,2,3]. However, applications should target complex but parsimonious models which make realistic assumptions and let the observed data drive estimations.

There are two common approaches to epidemiological modeling: phenomenological models and mechanistic models (e.g., compartmental models). On the one hand, phenomenological models use an empirical approach based on growth curve fitting (e.g., by nonlinear least squares [4] or by maximum likelihood [5]) to describe the temporal progression of case counts (e.g., daily confirmed positive cases). In this regard, the logistic bell curve has been widely used for various epidemic data, but it lacks flexibility for epidemics whose data exibits asymmetry or varying growth patterns [4,6,7]. With a view to allow flexibility, Tovissodé et al. [5] considered the generic growth curve of Turner et al. [8] with application to COVID-19 data. This approach concedes the simple logistic curve when it is supported by the observed data, but offers the possibility to fit various flexible growth models such as the generalized logistic model [9,10], the hyperlogistic model [8,11], the hyper-Gompertz [8] and the Gompertz curves [12,13]. However, to be realistic, models for epidemic data should be able to account for the potential effect of containment measures when implemented after an epidemic outbreak. In a target population undergoing an epidemic wave, the number of infective individuals may be assumed to follow an exponential growth in the early epidemic phase where no containment measures were implemented or the implemented measures were not yet effective [14]. In this case, the variation of the number of infective individuals is expected to shift to a sub-exponential growth resulting from negative feedbacks due to a decrease in the probability that an infectious individual meets a susceptible individual [6] or effects of the containment measures, if any. The major advantage of the phenomenological modeling approach is its simplicity while allowing the estimation of various quantities of interest to understand an epidemic, e.g., the “epidemic latency period” defined as the delays between the appearance of the first infectious case in the population (“patient zero”) and the outbreak [14] and epidemic peak time and size, and the forecast of future incidence. The main limitation of phenomenological models is the inability to inform on the transmission process (new infections) and the removal processes (recovery and death) of an epidemic. As a result, phenomenological modeling lacks the ability to assess the effects of control interventions.

On the other hand, and contrary to phenomenological models, mechanistic models structure the population under study into different epidemiological states [4] and allow assessing the effects of control interventions on the population and disease dynamics. For instance, the effect of various control measures (e.g., contact limitation, detection and diagnosis) on COVID-19 transmission has been assessed using the Susceptibles–Exposed–Infectives–Recovered (SEIR) model and its variants [15,16,17]. However, because only a few epidemiological states can be observed, mechanistic models often face an identifiability issue in the estimation of model parameters [18,19,20]. In addition, there is generally no closed form solutions to the differential equations describing the considered epidemiological states. As a consequence, the estimation of compartmental models often relies on numerical approximations which make fitting procedures (e.g., nonlinear least squares or Bayesian estimation) computationally intensive and may introduce high-order errors in both estimates and forecasts [21]. Moreover, some quantities of high interest to understand epidemic outbreaks, which are readily available from a growth model including the epidemic latency period, are hard to derive under compartmental models.

This study proposes a hybrid framework to combine the advantages of phenomenological and mechanistic models while circumventing some of the limits of the two approaches. We focus on epidemic waves managed with at least an isolation measure for all identified infectives, as for the COVID-19 pandemic in nearly all the world. The objective of this work is to provide a quantitative framework in which epidemiologists can identify, from a large family of models, the parsimonious model that explains patterns in an observed dataset, and then assess hypotheses on the potential course of related but unobservable processes of interest. Specifically, we modeled confirmed positive cases using a combination of the exponential growth curve for the initial epidemic phase and the generic growth curve [8] after this initial phase. This development allows the estimation of the duration of the exponential growth phase and the theoretical time and size of the peak of new positive cases. Secondly, we modeled removal (recovery and death) from identified positive cases as binary processes using two logistic regression models to monitor the evolution and the peak (time and size) of the actives among detected cases. Finally, to provide an overall view for a target epidemic, we integrated the growth curve and the logistic regression removal rates into a mechanistic SIQR model frame [22] in which the population is structured in Susceptibles, Infectives, Quarantined (identified actives cases) and Recovered individuals. The result is a mechanistic model in which the sizes of the different states (compartments) have closed form expressions. This allows inference on various epidemiological parameters such as the delay between the appearance of the first infectious case in the population (“patient zero”) and the outbreak (“epidemic latency period”), the reproduction number, the unobservable new infections per unit time as well as the proportion of the target population immunized against the pathogen of the target disease.

In addition to the estimates (with quantified uncertainty) for common epidemiological parameters, the proposed hybrid modeling framework extracts from the observed data and demographic rates, the evolution along the epidemic course of the key parameter to summarize the dynamics of an epidemic: the reproduction number. The changes in this parameter can thus be confronted to control measures promoted/enforced by public health authorities and governments. For illustrative purpose, we used the developed modeling framework: (i) to model COVID-19 case reporting data (daily PCR-confirmed positives, recoveries and deaths) from Western Africa (28 February to 31 August 2020); and (ii) to evaluate the transmission pattern of the disease in the region during the considered period. The results were used to discuss the effectiveness of some containment measures implemented by governments across the region.

2. The Hybrid Modeling Framework

In this section, we describe the three sub-models integrated into the proposed modeling framework, namely, the growth model, the logistic removal rates and the Susceptible–Infective–Quarantined–Recovered (SIQR) mechanistic model.

2.1. Mixture of Growth Models for Detected Cases

We assume that the cumulative number $C_t$ of reported cases, as a function of time t, has the form

$$\begin{array}{ccc}\hfill C_t& =& \left\{\begin{array}{cc}0& \mathrm{if}t\le 0\\ {\mathrm{e}}^{{\omega }_0(t-{\tau }_0)}& \mathrm{if}0<t\le t_e\\ \xi +{\phi }_t& \mathrm{if}t>t_e\end{array}\right.\hfill \end{array}$$

(1)

where $t_e>0$ is the duration from outbreak to the end of the exponential growth phase,

$$\begin{array}{c}\hfill {\phi }_t=\mathsf{\Omega }{(1+u_t)}^{-1/\nu }\end{array}$$

(2)

is the generic growth model [8] with $u_t={\left[1+\omega \nu \rho (t-\tau )\right]}^{-1/\rho }$, $\mathsf{\Omega }>0$ is a constant such that the ultimate epidemic size (detected) is $\xi +\mathsf{\Omega }$, $\omega >0$ is the “intrinsic” growth rate constant for the sub-exponential growth phase, $\nu >0$ is a growth acceleration parameter, $\rho$ ($-1<\rho <{\nu }^{-1}$) is a shape parameter controlling the skewness of the growth curve during the sub-exponential epidemic phase (see Appendix A.1 for restriction related details) and $\tau$ is a constant of integration determined by the initial conditions of the epidemic. The generic growth curve ${\phi }_t$ specified for $t>t_e$ encompasses many special or limiting cases including the Bertalanffy–Richards ($\rho \to 0$), hyper-Gompertz ($\nu \to 0$ while $\omega {\nu }^{1+\rho }\to \tilde{\omega }$ with $\tilde{\omega }$ constant), Gompertz ($\nu \to 0$, $\rho \to 0$ while $\omega \nu \to \tilde{\omega }$), hyper-logistic ($\nu =1$) and logistic ($\nu =1$ and $\rho \to 0$) growth models [8] (see Appendix A.1 for details). The parameter ${\omega }_0>0$ in (1) is the exponential growth rate for the early epidemic phase and ${\tau }_0\in \mathbb{R}$ determines the growth rate at $t=0$. The constants ${\omega }_0$ and ${\tau }_0$ are set such that the first derivative ${\dot{C}}_t$ and the second derivative ${\ddot{C}}_t$ of $C_t$ with respect to t are smooth at $t=t_e$ (i.e., at the end of the exponential growth phase). Specifically,

$$\begin{array}{ccc}\hfill {\omega }_0& =& {\ddot{\phi }}_e\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill {\tau }_0& =& t_e+\frac{log{\omega }_0-log{\dot{\phi }}_e}{{\omega }_0}\hfill \end{array}$$

(4)

where ${\dot{\phi }}_e={\dot{\phi }}_{t_e}$ and ${\ddot{\phi }}_e={\ddot{\phi }}_{t_e}$; ${\dot{\phi }}_t$ and ${\ddot{\phi }}_t$ are, respectively, the first and second derivatives of ${\phi }_t$ (see Appendix A.1 for details); and (4) follows from setting ${\omega }_0{\mathrm{e}}^{{\omega }_0(t_e-{\tau }_0)}={\dot{\phi }}_e$. Furthermore, the real constant $\xi$ in (1) ensures that $C_t$ does not jump at $t=t_e$. In other words, $\xi$ is given by $\xi ={\mathrm{e}}^{{\omega }_0(t_e-{\tau }_0)}-{\phi }_e$ (with ${\phi }_e={\phi }_{t_e}$) which by (4) simplifies to

$$\begin{array}{c}\hfill \xi =\frac{\dot{{\phi }_e}}{{\omega }_0}-{\phi }_e.\end{array}$$

(5)

In (1), the time (in e.g., days, weeks or months) of the first identified cases corresponds to $t=1$. In other words, to match (1) to the observed data, $C_1$ is identified to the number of cases reported in the time interval $(0,1]$, $C_2$ is the number of cases reported in the time interval $(0,2]$, etc. If $\mathsf{\Omega }\to \infty$ and $\nu \rho \to 0$, the curve $C_t$ converges to an exponential growth curve with rate ${\omega }_0$. However, this scenario can be ruled out since the size of any target population is finite and so is $\mathsf{\Omega }$. In practice, the exponential growth is prevented by negative feedbacks which decrease the probability that an infectious individual and a susceptible individual meet and have an adequate contact (i.e., contact sufficient for transmission). For instance, the growth of the infectives is naturally continuously lowered by the increasing fraction of the population constituted by individuals who recovered and become less susceptible (temporarily or permanently immune) to the infection [6]. To prevent the exponential growth of the infectives, control measures such as quarantining and lockdown reduce the probability of contact between susceptible and infectious individuals, whereas some other measures such as social distancing and wearing a face mask reduce the likelihood of transmission whenever contacts happen.

The specification of the growth model in (1) to an epidemic thus implies that the growth rate ${\dot{C}}_t$, i.e., the number of new cases reported per unit time given by

$$\begin{array}{c}\hfill {\dot{C}}_t=\left\{\begin{array}{cc}{\omega }_0{\mathrm{e}}^{{\omega }_0(t-{\tau }_0)}& \mathrm{if}0\le t\le t_e\\ {\dot{\phi }}_t& \mathrm{if}t>t_e\end{array}\right.\end{array}$$

(6)

with ${\dot{\phi }}_t$ defined in Appendix A.1, will peak and then fall toward zero case per unit time. The peak occurs at a time $t_p>t_e$ when the growth acceleration ${\ddot{C}}_t$ given by,

$$\begin{array}{c}\hfill {\ddot{C}}_t=\left\{\begin{array}{cc}{\omega }_0^2{\mathrm{e}}^{{\omega }_0(t-{\tau }_0)}& \mathrm{if}0\le t\le t_e\\ {\ddot{\phi }}_t& \mathrm{if}t>t_e\end{array}\right.\end{array}$$

(7)

with ${\ddot{\phi }}_t$ defined in Appendix A.1, vanishes. The expressions of the time ($t_p$) and the size (${\dot{C}}_p$) of the peak are available in Appendix A.2 for the general situation ($\nu \ne 0$ and $\rho \ne 0$), as well as for limiting cases.

The number of detected cases $C_t$ is the basic data reported during an epidemic. Once this has been modeled, various epidemic related quantities can be inferred upon introduction of disease related parameters (e.g., detection of infectives, recoveries and deaths) and demographic parameters (e.g., natural mortality, births and immigration).

2.2. Infectives, Epidemic Latency Period and Active Cases

Since only a fraction of infectives are identified at a time t, the number $I_t$ of infective individuals in a target population is obtained using (6) as $I_t={\delta }^{-1}{\dot{C}}_t$ [5], which reads

$$\begin{array}{ccc}\hfill I_t& =& \left\{\begin{array}{cc}{I}_0{\mathrm{e}}^{{\omega }_{0}t}& \mathrm{if}t\le t_e\\ {\delta }^{-1}{\dot{\phi }}_t& \mathrm{if}t>t_e\end{array}\right.\hfill \end{array}$$

(8)

where $I_0={\delta }^{-1}{\omega }_0{\mathrm{e}}^{-{\omega }_0{\tau }_0}$ is the number of infectives at the outbreak ($t=0$) and $\delta \in (0,1]$ is the detection rate assumed constant along the epidemic course (after the outbreak). Note that the number of infectives before the outbreak ($t<0$) is obtained by back extrapolation as $I_t=I_0{\mathrm{e}}^{{\omega }_{0}t}$, i.e., considering an exponential growth before the outbreak [14].

We refer to the time from the appearance of the first infectious case in the population (“patient zero”) to the outbreak as the “epidemic latency period”. An estimate of the duration $t_o$ of this period is obtained by setting $I_t=1$ [14]. By (8), the duration of the epidemic latency period is estimated by $t_o={\omega }_0^{-1}log{I}_0$, which on using (4) simplifies to

$$\begin{array}{c}\hfill t_o=\frac{log{\dot{\phi }}_e-log\delta }{{\omega }_0}-t_e.\end{array}$$

(9)

The number of detected and active cases, i.e., individuals tested positive and in isolation at a hospital or at home at time t, is denoted $Q_t$ following Hethcote et al. [22] for “Quarantined” state, although we refer to $Q_t$ as “Actives”. Given the detected cases $C_t$ in (1), $Q_t$ satisfies

$$\begin{array}{c}\hfill {\dot{Q}}_t={\dot{C}}_t-({\alpha }_t+{ϵ}_t)Q_t\end{array}$$

(10)

where ${\alpha }_t$ is the recovery rate and ${ϵ}_t$ is the death rate (natural and disease-related mortality) of actives. Indeed, following Tovissodé et al. [5], we allow the removal rates ${\alpha }_t$ and ${ϵ}_t$ from $Q_t$ to be time varying. This is appropriate when recovery and death data are available in addition to the reported positive cases per unit time. The two rates have here the logistic forms

$$\begin{array}{ccc}\hfill {\alpha }_t& =& {\left[1+{\mathrm{e}}^{-({\kappa }_0+\kappa t)}\right]}^{-1}\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill {ϵ}_t& =& {\left[1+{\mathrm{e}}^{-({\lambda }_0+\lambda t)}\right]}^{-1}.\hfill \end{array}$$

(12)

The number of active cases is then given by (see Appendix B for details)

$$\begin{array}{ccc}\hfill Q_t& =& \left\{\begin{array}{cc}\left[Q_0{F}_0+{\omega }_0{\int }_0^t{\mathrm{e}}^{{\omega }_0(r-{\tau }_0)}{F}_{r}dr\right]F_t^{-1}& \mathrm{if}0<t\le t_e\\ \left[Q_e{F}_{t_e}+{\int }_{t_e}^t{\dot{\phi }}_r{F}_{r}dr\right]F_t^{-1}& \mathrm{if}t>t_e\end{array}\right.\hfill \end{array}$$

(13)

where $Q_0$ is available from Equation (A3) and represents the number of persistent cases from previous epidemic waves (isolated actives) at the outbreak of the target epidemic wave (e.g., $Q_0=0$ for a new disease-related epidemic) and $F_t$ is defined as

$$\begin{array}{c}\hfill F_t=\left\{\begin{array}{cc}{\mathrm{e}}^{({\alpha }_0+{ϵ}_0)t}& \mathrm{if}\kappa =0\mathrm{and}\lambda =0\\ {\mathrm{e}}^{{\alpha }_{0}t}{\left(1+{\mathrm{e}}^{{\lambda }_0+\lambda t}\right)}^{1/\lambda }& \mathrm{if}\kappa =0\mathrm{and}\lambda \ne 0\\ {\left(1+{\mathrm{e}}^{{\kappa }_0+\kappa t}\right)}^{1/\kappa }{\mathrm{e}}^{{ϵ}_{0}t}& \mathrm{if}\kappa \ne 0\mathrm{and}\lambda =0\\ {\left(1+{\mathrm{e}}^{{\kappa }_0+\kappa t}\right)}^{1/\kappa }{\left(1+{\mathrm{e}}^{{\lambda }_0+\lambda t}\right)}^{1/\lambda }& \mathrm{if}\kappa \ne 0\mathrm{and}\lambda \ne 0\end{array}\right..\end{array}$$

(14)

2.3. Overall Epidemic Dynamics

The dynamics of an epidemic, as expressed by the variations of the infectives $I_t$, is determined by the combination of the transmission rate (new infections) and the average residence time, i.e., the average duration from infection to isolation, recovery or death. The core parameter to summarize these dynamics is, at moment t, the reproduction number denoted ${\mathcal{R}}_t$, which is indeed crucial for quantifying the intensity of control measures required to control an epidemic [7].

The reproduction number is defined as the average number of secondary cases generated by a primary case. With a view to derive ${\mathcal{R}}_t$ under the growth model in (1), we first consider an overall picture of the target population in order to enlighten the sources (transmission and removal) of the variations of $I_t$ as given in (8).

2.3.1. The SIQR Model

Following the authors of [5,14], we consider the Susceptible–Infectious–Quarantined–Recovered (SIQR) model of Hethcote et al. [22] to obtain a picture of the different states of individuals in a target population. We use the “quarantine-adjusted incidence” version [22] of this model since the underlying transmission mechanism explicitly recognizes the isolation of detected cases. In this framework, letting $N_t$ denote, at time t, the size of the target population (assumed finite but large), $N_t$ satisfies

$$\begin{array}{c}\hfill N_t=S_t+I_t+Q_t+R_t\end{array}$$

(15)

where $S_t$ is the size of the class of susceptible individuals, $I_t$ is the class of infectives, $Q_t$ is the size of the class of detected active cases and $R_t$ is the size of the class of individuals who recovered (both detected and not detected). We assume that the infection has zero latent period (susceptible individuals become infectious as soon as they become infected). The individuals in the classes R are assumed permanently immune within the period of time considered. It is also assumed that known active cases (in the class Q) do not mix with other classes and do not infect the susceptibles (i.e., the transmission rate from Q-class individuals is considered negligible). The corresponding SIQR model is described by the following set of nonlinear differential equations [22]

$$\begin{array}{ccc}\hfill {\dot{S}}_t& =& \eta -{\beta }_t(S_t+R_t)I_t/(N_t-Q_t)-\mu S_t\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill {\dot{I}}_t& =& \left[{\beta }_t(S_t+R_t)/(N_t-Q_t)-(\gamma +{\delta }_t+\pi )\right]I_t\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill {\dot{Q}}_t& =& {\delta }_t{I}_t-({\alpha }_t+{ϵ}_t)Q_t\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill {\dot{R}}_t& =& \gamma I_t+{\alpha }_t{Q}_t-\mu R_t\hfill \end{array}$$

(19)

where $\eta$ is the recruitment rate of susceptibles (births and immigration); ${\beta }_t$ is the total number of adequate contacts (i.e., contacts sufficient for transmission) per unit time; $\mu$ is the per capita natural mortality rate; ${\alpha }_t$ and $\gamma$ are the recovery rates from actives $Q_t$ and infectives $I_t$ respectively; ${ϵ}_t$ and $\pi$ are the death rates (natural and disease-related) for actives $Q_t$ and infectives $I_t$ respectively; and ${\delta }_t$ is the detection rate which is null (${\delta }_t=0$) for $t<0$ and equals ${\delta }_t=\delta$ for $t\ge 0$. Note that (18) is the same as (10) for $t\ge 0$. Unlike in [22], we allow the transmission rate ${\beta }_t$ to be time varying as a consequence of the form of the number of infectives $I_t$ already available in (8). The transfer diagram for this SIQR model is shown in Figure 1.

The system (16)–(19) always has the disease-free equilibrium $P_0=(S=\eta /\mu ,I=0,Q=0,R=0)$, i.e., in the absence of the disease, the population size $N_t$ approaches the carrying capacity $N^{*}=\eta /\mu$. Further discussion of the equilibria of the system are given in Appendix C.1. The availability of the number of infectives in Equation (8) makes it possible to solve the system (16)–(19). Indeed, from (17), the transmission rate, i.e., the number of adequate contacts per unit time (for $I_t>0$) is given by

$$\begin{array}{c}\hfill {\beta }_t=\left(\gamma +{\delta }_t+\pi +\frac{{\dot{I}}_t}{I_t}\right)\left(1+\frac{I_t}{S_t+R_t}\right).\end{array}$$

(20)

From (20), and using the same approach considered to find the number $Q_t$ of active cases in Equation (13) from the number $I_t$ of infectives in Equation (8), the expressions of the number $S_t$ of susceptibles, the number $R_t$ of recovered individuals and the total number of persons infected during an epidemic wave can be obtained (see Appendix C.2 for details).

2.3.2. The Effective Reproduction Number

From the definition of the effective reproduction number as the average number of secondary cases generated by a primary case, the threshold ${\mathcal{R}}_t$ corresponds to the product of the transmission rate ${\beta }_t$ and the average residence time $1/(\gamma +{\delta }_t+\pi )$ in the class of infectives, i.e.,

$$\begin{array}{c}\hfill {\mathcal{R}}_t={\beta }_t/(\gamma +{\delta }_t+\pi ).\end{array}$$

This effective reproduction number is sometimes referred to as a “quarantine” reproduction number [22] or simply a “control” reproduction number to acknowledge the influence of isolation of identified infectives, and other control measures, if any [15]. The basic reproduction number defined as the average number of secondary infections produced when one primary infectious individual enters a completely susceptible population ($S_o=N_o-1$, $I_o=1$, $Q_o=0$, $R_o=0$), is here given by ${\mathcal{R}}_o=\left(1+\frac{{\omega }_0}{\gamma +\pi }\right)N_o/(N_o-1)$. This expression is simplified, assuming $N_o/(N_o-1)=1$ for the sake of beauty [23] and mostly because $N_o$ is large (recall this is a model assumption), as

$$\begin{array}{c}\hfill {\mathcal{R}}_o=1+\frac{{\omega }_0}{\gamma +\pi }.\end{array}$$

(21)

During the epidemic latency period ($t_o<t<0$) where the growth is exponential (${\dot{I}}_t/I_t={\omega }_0$) and the detection rate is ${\delta }_t=0$, the time-varying reproduction number is given by

$$\begin{array}{c}\hfill {\mathcal{R}}_t=\left(1+\frac{{\omega }_0}{\gamma +\pi }\right)\left(1+\frac{I_t}{S_t+R_t}\right)\mathrm{for}-t_o\le t<0.\end{array}$$

(22)

From the outbreak, the time-varying effective reproduction number during the remaining of the exponential phase has the same form

$$\begin{array}{c}\hfill {\mathcal{R}}_t=\left(1+\frac{{\omega }_0}{\gamma +\delta +\pi }\right)\left(1+\frac{I_t}{S_t+R_t}\right)\mathrm{for}0\le t\le t_e.\end{array}$$

(23)

It appears from (22) and (23) that ${\mathcal{R}}_t>1$ during the whole exponential growth phase as expected. During the sub-exponential growth phase, the time-varying effective reproduction number is given by

$$\begin{array}{c}\hfill {\mathcal{R}}_t=\left(1+\frac{z_t}{\gamma +\delta +\pi }\right)\left(1+\frac{I_t}{S_t+R_t}\right)\mathrm{for}t>t_e\end{array}$$

(24)

where $z_t={\ddot{\phi }}_t/{\dot{\phi }}_t$ (see Appendix A.1).

2.3.3. Epidemic Peak

The peak of new infections occurs when the second derivative of the total number of infected persons (since the beginning of the epidemic) vanishes. This peak time denoted $t_{new}$ satisfies $t_{new}>t_e$ and is the solution of (see details in Appendix C.3)

$$\begin{array}{c}\hfill (\gamma +\delta +\pi ){\ddot{\phi }}_t+{\stackrel{⃛}{\phi }}_t=0\end{array}$$

(25)

which can be solved for t using a numerical root finding routine such as the R [24] function uniroot or the Matlab [25] function fzero. Afterwards, the peak size ${\dot{T}}_{new}$ (the maximum number of new infections per unit time) is obtained by inserting $t_{new}$ in (A14).

2.4. Long-Term Epidemic Dynamics

The specification of the growth model in (1) to an epidemic implicitly assumes that the number of infectives in (8) peaks at time $t_p$ and then approaches zero. The decay of the infectives after the peak can happen at various rates, depending on the growth pattern (determined by contacts between the infectives and the susceptibles or intermediate hosts), the response of the infected individual’s organism (natural or induced with medicine or a vaccine) to the disease (recovery and death process) and the testing efforts (detection followed by isolation). There are actually two alternative paths from a disease-related state (i.e., $I_t>0$) toward the unique (disease-free) equilibrium $P_0$: transmissions either stop (${\mathcal{R}}_t$ reaches zero) or continue fro a long time at a rate which cannot sustain an epidemic ($0<{\mathcal{R}}_t\le 1$). These two scenarios are discussed further in Appendix C.4.

2.5. Statistical Model and Inference

To allow likelihood inference in the growth models in (1) using observed epidemiological data, we follow Tovissodé et al. [5] and assign to new reported cases $Y_t$ ($t=1,2,\cdots ,n$) a log-normal distribution with probability density function (pdf)

$$\begin{array}{c}\hfill f_Y\left(Y_t\right|\mathbf{\theta })=\frac{1}{\sigma (Y_t+1)\sqrt{2\pi }}exp\left(-\frac{1}{2}{\left[\frac{log(Y_t+1)-log({\dot{C}}_t+1)}{\sigma }+\frac{\sigma }{2}\right]}^2\right)\end{array}$$

(26)

where $\sigma >0$ is a dispersion parameter (standard deviation at logarithmic scale). This specification yields the mean $E\left[Y_t\right]={\dot{C}}_t$ and the variance $Var\left[Y_t\right]={\left({\dot{C}}_t+1\right)}^2\left({\mathrm{e}}^{{\sigma }^2}-1\right)$ while allowing null values of $Y_t$. In addition, the numbers of new recoveries $G_t$ and new deaths $M_t$ from known active cases $Q_t$ ($t=1,2,\cdots ,n$) are modeled using logistic regression models with probability mass functions (pmf)

$$\begin{array}{ccc}\hfill f_G\left(G_t\right|\mathbf{\theta },Q_{t-1},Y_t)& =& \left(\genfrac{}{}{0pt}{}{Q_{t-1}+Y_t}{G_t}\right){\alpha }_t^{G_t}{(1-{\alpha }_t)}^{Q_{t-1}+Y_t-G_t}\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill f_M\left(M_t\right|\mathbf{\theta },Q_{t-1},Y_t)& =& \left(\genfrac{}{}{0pt}{}{Q_{t-1}+Y_t}{M_t}\right){ϵ}_t^{M_t}{(1-{ϵ}_t)}^{Q_{t-1}+Y_t-M_t}\hfill \end{array}$$

(28)

where ${\alpha }_t={\left[1+{\mathrm{e}}^{{\kappa }_0+\kappa t}\right]}^{-1}$ and ${ϵ}_t={\left[1+{\mathrm{e}}^{{\lambda }_0+\lambda t}\right]}^{-1}$. The parameter vector indexing the pdf in (26) and the conditional pmf in (27) and (28) is $\mathbf{\theta }={(\mathsf{\Omega },\omega ,\nu ,\rho ,\tau ,t_e,\sigma ,{\kappa }_0,\kappa ,{\lambda }_0,\lambda )}^{\top }$ when the generic growth curve is considered for the sub-exponential growth phase. If a special case of the generic growth curve is desired, the corresponding restricted parameters must be withdrawn from $\mathbf{\theta }$. For instance, the use of a hyper-logistic growth curve ($\nu =1$) implies $\mathbf{\theta }={(\mathsf{\Omega },\omega ,\rho ,\tau ,t_e,\sigma ,{\kappa }_0,\kappa ,{\lambda }_0,\lambda )}^{\top }$. Given $Q_0$, the conditional log-likelihood of an observed series $\{Y_t,G_t,M_t\}$ with $t=1,2,\cdots ,n$, as a function of the parameter $\mathbf{\theta }$ is

$$\begin{array}{ccc}\hfill \ell \left(\mathbf{\theta }\right)& =& {\ell }_Y\left(\mathbf{\theta }\right)+{\ell }_G\left(\mathbf{\theta }\right)+{\ell }_M\left(\mathbf{\theta }\right)\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill \mathrm{where} {\ell }_Y\left(\mathbf{\theta }\right)& =& \sum _{t=1}^{n}log{f}_Y\left(Y_t\right|\mathbf{\theta })\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill {\ell }_G\left(\mathbf{\theta }\right)& =& \sum _{t=1}^{n}log{f}_G\left(G_t\right|\mathbf{\theta },Q_{t-1},Y_t)\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill {\ell }_M\left(\mathbf{\theta }\right)& =& \sum _{t=1}^{n}log{f}_M\left(M_t\right|\mathbf{\theta },Q_{t-1},Y_t).\hfill \end{array}$$

(32)

The conditional maximum likelihood estimate $\widehat{\mathbf{\theta }}$ of $\mathbf{\theta }$ can be obtained using an optimization algorithm to maximize the log-likelihood function ℓ. Note that the three components of $\ell \left(\mathbf{\theta }\right)$ are separable and can be maximized independently. In other words, the parameter vector $\mathbf{\theta }$ has the partition $\mathbf{\theta }={({\mathbf{\theta }}_Y^{\top },{\mathbf{\theta }}_G^{\top },{\mathbf{\theta }}_M^{\top })}^{\top }$ and the maximum likelihood estimates of the components ${\mathbf{\theta }}_Y={(\mathsf{\Omega },\omega ,\nu ,\rho ,{\tau }_0,t_e,\sigma )}^{\top }$, ${\mathbf{\theta }}_G={({\kappa }_0,\kappa )}^{\top }$ and ${\mathbf{\theta }}_M={({\lambda }_0,\lambda )}^{\top }$ can be obtained by maximizing ${\ell }_Y$, ${\ell }_G$ and ${\ell }_M$ respectively.

Since both the binomial and the log-normal distributions belong to the exponential family, we consider the common deviance statistic used in Generalized Linear Models [26] for checking the goodness-of-fit of the log-normal model associated to $Y_t$ and the binomial models associated to $G_t$ and $M_t$. For the selection of the parsimonious model agreeing with the observed data, we consider the likelihood ratio statistic [27]. Further details on the deviance statistic and the likelihood ratio test are given in Appendix D.

3. Application to COVID-19 Data of Western Africa

3.1. Context and Objectives

The Western African region has 16 countries (Benin, Burkina-Faso, Cape Verde, Côte d’Ivoire, Gambia, Ghana, Guinea, Guinea-Bissau, Liberia, Mali, Mauritania, Niger, Nigeria, Senegal, Sierra Leone and Togo), covering 6,140,178 km² with a population of about 402,555,230 inhabitants [28] (

Table 1. Population size [28] and cumulative PCR-confirmed COVID-19 cases, deaths and recoveries in West Africa (28 February to 31 August 2020) [29].

Country	Population Size	Total Confirmed	Recoveries	Deaths
Nigeria	206,522,290	54,008	41,638	1013
Ghana	31,072,945	44,298	42,963	276
Côte d’Ivoire	26,428,999	18,067	16,699	117
Niger	24,269,389	1176	1088	69
Burkina-Faso	20,946,992	1368	1058	55
Mali	20,294,900	2776	2169	126
Senegal	16,776,618	13,611	9439	284
Guinea	13,160,021	9409	8447	59
Benin	12,123,200	2145	1738	40
Togo	8,293,924	1400	1005	28
Sierra Leone	7,989,949	2022	1594	70
Liberia	5,066,990	1304	872	82
Mauritania	4,659,052	7048	6464	159
Gambia	2,421,823	2963	1032	96
Guinea-Bissau	1,971,640	2205	1127	34
Cape Verde	556,498	3884	2916	40
West Africa	402,555,230	167,684	140,249	2548

).

The first COVID-19 patient was formally identified in Western Africa in late (27) February 2020. We considered COVID-19 daily infection (PCR-confirmed cases on the day of reporting), recovery and death data, from 28 February to 31 August 2020, obtained from the Global Rise of Education Platform [29]. This period roughly corresponds to the first wave of the COVID-19 pandemic in the region [30]. We concentrated on these six months of data since the proposed modeling framework has been designed for a single epidemic wave. As of 31 August 2020, the region had 167,684 confirmed cases, among which $83.64\%$ recovered and $1.52\%$ died (Table 1). Although the region is heterogeneous, we treated it as if it were homogeneous. Indeed, it must be kept in mind that the reported COVID-19 cases occurred in small clusters concentrated in the main cities of each country. Hence, the sparsity of the data for the whole region actually reflect data sparsity at national and city levels.

The purpose of this analysis is to demonstrate, by example, the use of the proposed modeling framework. The specific aims are: (i) to model COVID-19 case reporting data (daily PCR-confirmed positives, recoveries and deaths) from Western Africa (28 February to 31 August 2020); and (ii) to evaluate the transmission pattern of the disease. Most West African governments have planned and subsequently implemented several control measures, either before or overlapping with the time of diagnosis of the first national cases [31]. The main sequence of public health and movement restriction measures taken by West African governments during the considered period includes personal hygiene and social distancing recommendations and isolation/lockdown (

Table 2. Main sequence of public health and movement restriction measures taken by West African governments during the first phase of the COVID-19 pandemic (until 31 August 2020).

Main Interventions	First Introduction (Country)	Implementation by the Last Country
State of health emergency and social distancing	22 March 2020 (Ghana)	30 March 2020 (Sierra Leone)
Setting up test sites and measures to quarantine suspected cases and isolate positive cases	25 February 2020 (Nigeria)	Early March 2020
Partial lockdown	18 March 2020 (Benin)	Late March
Curfew	20 March 2020 (Burkina Faso)	Not all countries
Reduced mobility and prohibition of social gatherings	15 March 2020 (Ghana)	Late March 2020
Land borders closure	20 March 2020 (Côte d’Ivoire)	30 March 2020 (Sierra Leone)
Wearing face mask in public mandatory	8 April 2020 (Benin)	14 May 2020 (Mauritania)
Systematic testing of target groups	22 March 2020 (Benin)	Late March to early April

Sources: https://hsfnotes.com/africa/2020/05/22/covid-19-initial-responses-of-certain-african-countries#page=1 (accessed on 4 April 2021) and [31].

). The adoption of these containment measures followed a sustained increment during late March 2020. The modeling results are used to discuss the effectiveness of the containment measures and the implications for the control of the further spread of COVID-19 in West African countries.

3.2. Data Analysis

All computations and statistical analyses were performed in R software [24]. The significance level of statistical tests was set to 5%.

3.2.1. Model Fitting

We fitted the generic growth curve to the daily new infections $Y_t$. We used the optim routine of R software to maximize the log-likelihood (30). We also fitted three of its special cases (Bertalanffy–Richards, hyper-logistic and hyper-Gompertz), which were compared to the generic model fit using likelihood ratio tests. Instead of directly maximizing the log-likelihoods (31) for ${\widehat{\mathbf{\theta }}}_G$ and (32) for ${\widehat{\mathbf{\theta }}}_M$ with the optim routine, we used the glm routine of R with the family specification “family = binomial(logit)”. Since COVID-19 was a new disease in 2020, we considered the number of known active cases $Q_0=0$ at $t=0$ in (27) and (28). We plotted the daily new positives, recoveries, deaths and actives to provide graphical insights in the fitted models.

3.2.2. Overall Epidemic Dynamics

We analyzed the overall dynamics of the COVID-19 epidemic in West Africa using the mechanistic SIQR model described in Section 2.3. The rate parameters $\delta$ (detection rate), $\gamma$ and $\pi$ (recovery and death rates in infected but non-detected individuals) cannot be estimated using only the available data sequence $\{Y_t,G_t,M_t\}$ (daily new positives, recoveries and deaths) without additional assumptions on their relationships with the rate parameters for detected cases (${\alpha }_t$ and ${ϵ}_t$). We obtained from the literature $\delta =0.009$ [30] and $\gamma +\pi =1/10$ [14,30] and assumed that the ratio of the daily recovery probability to the daily death probability in non detected infectives is equal to this ratio in the detected individuals at outbreak, i.e., before the implementation of treatments, if any. From $\gamma /\pi ={\alpha }_0/{ϵ}_0\approx 5.1495$, we obtained $\gamma =1/11.9419$ and $\pi =1/61.4953$.

Two demographic parameters are required in the SIQR model: the daily recruitment rate of susceptibles (through births and immigration) $\eta$ (individuals/day) and the per capita natural mortality rate $\mu$ (day⁻¹). Using the birth rate ${\rho }_b$ (total births and net immigrations in a period of length L divided by the average population size $\overline{N}$ during this period), the recruitment rate $\eta$ was estimated by

$$\begin{array}{c}\hfill \eta =\frac{r_b\overline{N}}{L}.\end{array}$$

(33)

Under “natural” (i.e., disease-free) conditions where $N_t=S_t$, the variation $\Delta N$ of the population size $N_t$ over a period of length L satisfies

$$\begin{array}{ccc}\hfill \Delta N& =& \left(\frac{\eta }{\mu }-N_i\right)\left(1-{\mathrm{e}}^{-\mu L}\right)\hfill \end{array}$$

(34)

where $N_i$ is the population of West Africa at the beginning of the period. The Equation (34) follows by (A5) with $I_0=0$. The variation $\Delta N$ of the population size is given by $\Delta N=r_b\overline{N}-r_d\overline{N}$, where $r_b\overline{N}$ represents the total recruitment during the period and $r_d\overline{N}$ represents the total number of deaths with ${\rho }_d$ the mortality rate (individuals/day). Consequently, $\mu$ can be obtained by solving (34) for $\mu$ using $\Delta N=(r_b-r_d)\overline{N}$.

We considered $L=365.25$ days, $\overline{N}=\mathrm{401,861,254}$, $N_i=\mathrm{397,429,929}$ [28]. Using the annual birth (32.816/1000) and death rates (7.952/1000) [32] and the net annual immigrations (−177,000 individuals) in West Africa [28], we obtained the rates $r_b=(32.816/1000)-(\mathrm{177,000}/\overline{N})=32.371/1000$ and $r_d=7.952/1000$. By (33) and (34), we then found and used for our analyses on West Africa, $\eta =\mathrm{35,615.35}$ individuals/day and $\mu =2.1745\times {10}^{-5}$ day⁻¹. We plotted the daily number of new infections, infectives and recovered individuals, as well as the reproduction number in the West African population.

3.2.3. Standard Error and Confidence Interval

Standard errors were obtained for quantities calculated using estimated model parameters by the delta method [33]. For a positive definite parameter or calculated quantity $\varphi$ in general, we first found the estimate $\widehat{\varphi }$ and its logarithmic scale-standard error ${\widehat{\sigma }}_{\varphi }$ by the delta method and computed its logarithmic scale-mean given by ${\widehat{\mu }}_{\varphi }=log\widehat{\varphi }-0.5{\widehat{\sigma }}_{\varphi }^2$. We then obtained the bounds of its shortest confidence interval as described by Dahiya and Guttman [34].

3.3. Results

3.3.1. Growth Curve for New Positives and Logistic Regressions for Removals

The results of the likelihood ratio tests comparing the generic growth model (1) against its closest special cases are presented in

Table 3. Likelihood ratio test results comparing the generic growth model [8] to three of its special cases.

Special Growth Model	Restriction	LRS	DF	p-Value
Bertalanffy–Richards	$\rho \to 0$	60.06	1	<0.001
Hyper-logistic	$\nu =1$	240.33	1	<0.001
Hyper-Gompertz	$\nu \to 0$, $\nu {\omega }^{(1+\rho )}\to \tilde{\omega }$	512.91	1	<0.001

Table notes: LRS, likelihood ratio statistic; DF, Degrees of freedom; $\tilde{\omega }$ stands for a positive constant; see Equation (1) for details on the parameters $\nu$, $\rho$ and $\omega$.

. The growth model involving the generic growth curve was retained. Indeed, the combination of an early exponential growth and the generic growth models was found to be the best growth model for the new positive cases in West Africa, as compared to the combinations of the exponential growth with the Bertalanffy–Richards, hyper-logistic and hyper-Gompertz growth models (Table 3; p-value < 0.001).

The deviance based ${\chi }^2$ test for overall goodness-of-fit (

Table 4. Deviance based goodness-of-fit test results for the combination of an early exponential growth curve with a generic growth curve (fitted to daily PCR-confirmed positives) and logistic regression models (fitted to daily numbers of recoveries and deaths) using West African COVID-19 data from 28 February to 31 August 2020.

Data	Goodness-of-Fit			Overall Significance
	DS	DF	p-Value	${\mathit{r}}_{\mathit{dev}}^2$ (%)	LRS	DF	p-Value
Reported cases	173.04	179	0.6115	95.26	3601.57	6	<0.001
Recoveries	45,028.51	184	<0.001	9.25	4861.55	1	<0.001
Deaths	499.13	184	<0.001	49.08	486.40	1	<0.001
Overall	45,700.68	175	<0.001	11.60	8949.52	8	<0.001

Table notes: DS, Deviance statistic; DF, Degrees of freedom; $r_{\mathit{dev}}^2$, adjusted-deviance reduction ratio; LRS, Likelihood Ratio Statistic.

) indicates a lack-of-fit (p-value < 0.001), with an overall adjusted-deviance reduction ratio of $r_{dev}^2=11.60\%$. Looking for the sub-models, we noticed that the estimated growth curve is significantly different from the corresponding null model fit (p-value < 0.001) and does not lack fit (p-value = 0.6115). Indeed, the adjusted-deviance reduction ratio is $r_{dev}^2=95.26\%$ (the adjusted-coefficient of determination is $r_a^2=99.96\%$). The overall lack of fit is due to the logistic regression fits for the daily recoveries ($r_{dev}^2=9.25\%$) and deaths ($r_{dev}^2=49.08\%$). We nevertheless kept these fits because there are significantly different from the corresponding null model fits (p-value < 0.001).

The maximum likelihood estimates of the generic growth model and logistic regression model parameters are presented in

Table 5. Estimate, standard error ($SE$), Wald test statistic (z-value), p-value ($\mathrm{P}(>|\mathrm{z}\left|\right)$) and 95% confidence interval ($C{I}_{95\%}$) for the parameters of the combination of an early exponential growth curve with a generic growth curve (fitted to daily PCR-confirmed positives) and logistic regression parameters (fitted to daily numbers of recoveries and deaths) using West African COVID-19 data from 28 February to 31 August 2020.

Parameter	Estimate	$\mathit{SE}$	z-Value *	$\mathbf{P}(>|\mathbf{z}\left|\right)$	${\mathit{CI}}_{95\%}$
$t_e$ ($day$)	29.4781	1.2368	80.1417	<0.001	[26.9413, 31.7865]
$\mathsf{\Omega }$ ($ind.$)	191,290.8	6444.5420	360.9696	<0.001	[178,756.4, 204,008.2]
$\omega$ (${day}^{-1}$)	0.0148	0.0007	−87.1715	<0.001	[0.0134, 0.0162]
$\nu$	3.7640	0.5280	9.3782	<0.001	[2.7685, 4.8240]
$\rho$	0.1202	0.0169	−15.1710	<0.001	[0.0884, 0.1541]
$\tau \left(day\right)$	171.3210	2.4252	70.6431	<0.001	[166.5678, 176.0742]
$\sigma$ ($logind.$)	0.3962	0.0201	−18.4774	<0.001	[0.3572, 0.4361]
${\kappa }_0$	−4.0609	0.0122	−333.6829	<0.001	[−4.0848, −4.0370]
$\kappa$	0.0059	0.0001	68.5372	<0.001	[0.0058, 0.0061]
${\lambda }_0$	−5.7136	0.0682	−83.7346	<0.001	[−5.8473, −5.5799]
$\lambda$	−0.0126	0.0006	−22.4195	<0.001	[−0.0137, −0.0115]
${\omega }_o$ (${day}^{-1}$)	0.1660	0.0011	−261.8024	<0.001	[0.1659, 0.1662]
${\tau }_0$ ($day$)	−7.2208	0.0226	−319.1971	<0.001	[−7.2651, −7.1764]
$\xi$ ($ind.$)	200.3128	2.7771	382.2758	<0.001	[194.8864, 205.7716]

Table notes: $ind.$, $individuals$; $t_e$ (day) is the duration of the exponential growth phase after the outbreak; $\mathsf{\Omega }$ and $\xi$ ($ind.$) determine the ultimate epidemic size (detected) as $\xi +\mathsf{\Omega }$; $\omega >0$ (day⁻¹) is the “intrinsic” growth rate constant for the sub-exponential growth phase; $\nu >0$ is a growth acceleration parameter, $\rho$ is a shape parameter controlling the skewness of the growth curve during the sub-exponential growth phase; $\tau$ (day) is a constant of integration determined by the initial conditions of the epidemic outbreak; $\sigma$ is the logarithmic-scale standard deviation of the log-normal distribution fitted to the daily new positive case reporting data; ${\kappa }_0$ and $\kappa$ are the logit-scale intercept and slope for the daily probability ${\alpha }_t$ that an active case recovers at time t (${\alpha }_t=1/(1+{\mathrm{e}}^{-({\kappa }_0+\kappa t)})$); ${\lambda }_0$ and $\lambda$ are the logit-scale intercept and slope for the daily probability ${ϵ}_t$ that an active case dies at time t (${ϵ}_t=1/(1+{\mathrm{e}}^{-({\lambda }_0+\lambda t)})$); ${\omega }_0$ (day⁻¹), ${\tau }_0$ (day) and $\xi$ (individuals) are not free parameters, but computed using Equations (4) and (5); ${\omega }_0$ is the growth rate during the exponential growth phase; ${\tau }_0$ and $\xi$ ensure that the daily number of positives ${\dot{C}}_t$ and the cumulative number of positives $C_t$ are smooth at $t=t_e$; * z-value was computed at logarithmic scale for positive definite parameters ($t_e$, $\mathsf{\Omega }$, $\omega$, $\nu$, $\rho$, $\sigma$ and ${\omega }_0$), so that a p-value $<0.05$ indicates significant difference from 1 at 5% level.

. The Wald test results (Table 5) agree with the likelihood ratio tests considered to select the growth model for the sub-exponential growth phase. Indeed, the 95% confidence bounds for the parameters $\nu$ ($CI\left(\nu \right)=[2.77,4.82]$) and $\rho$ ($CI\left(\rho \right)=[0.09,0.15]$) indicate that none of the Bertalanffy–Richards growth model ($\rho \to 0$), the hyper-logistic growth model ($\nu =1$), the logistic growth model ($\rho \to 0$, $\nu =1$), the hyper-Gompertz growth model ($\nu \to 0$, $\omega {\nu }^{1+\rho }\to \tilde{\omega }$) and the Gompertz growth model ($\rho \to 0$, $\nu \to 0$, $\omega \nu \to \tilde{\omega }$) are appropriate for this dataset.

The exponential growth phase lasted about one month (${\widehat{t}}_e=29.48$, $CI\left(t_e\right)=[26.94,31.79]$ days) after the outbreak (Table 5). The growth curve fitted to the cumulative positive cases is given by

$$\begin{array}{c}\hfill C_t=\left\{\begin{array}{c}{\mathrm{e}}^{0.1660\times (t+7.2208)}\mathrm{if}0<t\le 29.48\\ 200.3128+\frac{\mathrm{191,290.8}}{{\left\{1+{\left[1+0.0067\times (t-171.3210)\right]}^{8.3185}\right\}}^{0.2656}}\mathrm{if}t>29.48\end{array}\right.\end{array}$$

(35)

where t is the time (day) from the outbreak. Figure 2A shows the daily confirmed positive cases and the fitted growth curve based on a log-normal error structure. The observed peak of new positives happened 148 days after the outbreak (24 July 2020) and amounted to 2626 positive cases. However, the number of positive cases showed a high variability around this date (16–29/07/2020), with most daily records roughly ranging between 1600 and 2000 new positive cases (Figure 2A) around an average of 1803 cases (with standard error $SE=86.48$). The estimated peak time for the new positive cases was around 15 July 2020, i.e., about 139 days after the outbreak (

Table 6. Estimate, standard error ($SE$) and 95% confidence interval ($C{I}_{95\%}$) for some quantities using the West African COVID-19 data from 28 February to 31 August 2020.

Quantity	Observed Value	Estimate	$\mathit{SE}$	${\mathit{CI}}_{95\%}$
$t_o$ ($day$)	-	24.78	2.55	[19.91, 29.87]
$R_o$	-	2.66	0.11	[2.60, 2.69]
$I_0$ ($ind.$)	-	61.17	6.94	[47.98, 75.05]
$R_0$	-	2.52	0.12	[2.29, 2.76]
$t_p$ ($day$)	148	138.87	2.26	[134.45, 143.31]
${\dot{C}}_p$ ($ind.$)	2626	1804.90	83.40	[1643.19, 1969.86]
$t_{new}$ ($day$)	-	131.12	2.53	[126.18, 136.11]
${\dot{T}}_{max}$ ($ind.$)	-	22,352.97	1067.46	[20,284.04, 24,464.98]
$t_{Q_{max}}$ ($day$)	143	149.67	1.78	[146.18, 153.17]
$Q_{max}$ ($ind.$)	41,435	42,507.01	1449.81	[39,687.48, 45,368.24]
$R_{186}$ ($ind.$)	-	1,754,698.5	40,665.66	[1,675,407.60, 1,834,783.00]

Table notes: $t_o$, duration of the epidemic latency period; $R_o$, basic reproduction number; $I_0$, number of infectives at outbreak; $R_0$, reproduction number at outbreak; $t_p$, time of the peak of positive cases; ${\dot{C}}_p$, size of the peak of positive cases; $t_{new}$, time of the peak of new infections; ${\dot{T}}_{max}$, size of the peak of new infections; $t_{Q_{max}}$, time of the peak of active cases; $Q_{max}$, size of the peak of active cases; $R_{186}$, total number of recovered in the population at $t=186$ days (i.e., at the end of the studied period (31 August 2020)); - indicates not applicable.

), and the estimate of the peak size is about 1805 new positive cases ($CI\left({\dot{C}}_p\right)=[1643.19,1969.86]$). Assuming a log-normal distribution, the 95% prediction interval for the peak size is $PI\left({\dot{C}}_p\right)=[1368.93,2669.55]$ new positive cases, which includes the observed value. The 95% prediction interval for the peak time is $PI\left(t_p\right)=[126.59,151.65]$ days, which also includes the observed peak time.

Based on the logistic regression parameters shown in Table 5, the probabilities of removals from the actives (quarantined) are shown in Figure 3. The probabilities of recovery and death are ${\widehat{\alpha }}_0=0.0169$ and ${\widehat{ϵ}}_0=0.0033$, respectively, at outbreak ($t=0$). The recovery probability then improved, with an odd ratio (recover/not recover) increasing on average by $0.59\%$ ($CI\left(\kappa \right)=[0.58,0.61]$%) each day. The death probability on the contrary decreased, with an odd ratio (die/not die) decreasing on average by $1.26\%$ ($CI\left(\lambda \right)=[-1.37,-1.15]$%) each day.

Figure 2B,C shows the removals (daily recovery and death) and the fitted values based on the logistic regression models for removal probabilities. We noticed that the lack-of-fit (indicated by the residual deviance test) is due to the very large variability of the observed daily proportions of recoveries and deaths. However, despite the lack-of-fit in the logistic regression fits, the use of the related recovery and death probabilities (${\alpha }_t$ and ${ϵ}_t$) along with the fitted growth curve (${\dot{C}}_t$), resulted in fitted active cases ($Q_t$) agreeing to a large extent with the observed daily actives (Figure 2D), with an adjusted-coefficient of determination of $97.08\%$. The peak of known active cases ($Q_t$) was on 19 July 2020 and amounted to 41,435 actives. The fitted peak is about 42,507 actives around 26 July 2020 (Table 6). The 95% prediction interval is $PI\left(Q_{max}\right)=[\mathrm{34,807.25},\mathrm{50,893.54}]$ actives for the maximum of active cases and $PI\left(t_{Q_{max}}\right)=[139.92,159.71]$ days for the peak time $t_{Q_{max}}$ (16 July to 5 August 2020).

3.3.2. Overall Epidemic Dynamics

The estimate of the duration of the epidemic latency period (delay between the arrival of the first infectious individual and outbreak) is about 25 days ($CI\left(t_o\right)=[19.91,29.87]$ days; see Table 6). Accordingly, the first imported COVID-19 case(s) in West Africa likely entered the region during the last week of January and the first week of February (28 January–7 February) 2020. The estimate of the basic reproduction number is ${\widehat{\mathcal{R}}}_o=2.66$ ($CI\left({\mathcal{R}}_o\right)=[2.60,2.69]$). At outbreak, the number of infectives in the region is estimated at about 61 ($CI\left(I_0\right)=[47.98,75.05]$) infectives. The estimate of the control reproduction number during the exponential growth phase after the outbreak is ${\widehat{\mathcal{R}}}_0=2.52$ ($CI\left({\mathcal{R}}_0\right)=[2.29,2.76]$).

Figure 4 shows the curves of the daily number of new infections (${\dot{T}}_t$), the daily number of infectives ($I_t$) and the immune fraction of the population ($R_t=K_t+U_t$). As expected, the peak in new infections occurred before the peak in detected infected individuals (observed 143 days after the outbreak). Indeed, the number of new infections peaked about 131 days after the outbreak ($CI\left(t_{new}\right)=[126.18,136.11]$ days), i.e., around 7 (2–12) July 2020, to about 22,353 ($CI\left({\dot{T}}_{max}\right)=[\mathrm{20,284.04},\mathrm{24,464.98}]$) new infections. As of 31 August 2020, the number of known recoveries in the West African region was 140,249. The number of both known and unknown recovered people at this date is estimated at about 1,754,699 individuals ($CI\left(R_{186}\right)=[\mathrm{1,675,407.60},\mathrm{1,834,783.00}]$), i.e., about 0.44% of the population in the region.

The time-varying effective reproduction number is shown in Figure 5. It appears that the effective reproduction number first decreased during the sub-exponential growth phase (from 2.52 on 27 February 2020), reaching 1 on 15 July and 0.66 on 31 August 2020. The effective reproduction number attained a minimum value of 0.61 on 29 September 2020 and then increased with a dynamics indicating ${\mathcal{R}}_{\infty }=1$.

4. Discussion

The importance of mathematical models in understanding and predicting the course of an epidemic outbreak and in assessing the impacts of public health control measures has been well documented in the current context of the COVID-19 pandemic [15,35,36,37]. Whereas phenomenological modeling is limited in the scope of inference, compartmental modeling faces identifiability issues and is usually computationally intensive [38]. This study proposes a hybrid modeling framework which combines phenomenological and mechanistic modeling approaches to assess the dynamics of epidemic outbreaks while circumventing some of the limitations of each approach. We illustrate our description of the different epidemiological aspects that the hybrid modeling framework deals with using COVID-19 data from West Africa (28 February to 31 August 2020). It is worth noting that the heterogeneity of the West African region in terms of testing and reporting policies, especially for the first epidemic wave, is an important limitation for this application. This is systematically true for any regional assessment of the pandemic [15]. Our analysis aims to provide an overall view of the dynamics of the pandemic in the West Africa. However, the analysis of the data from each country may be conducted to obtain finer country-specific results (for some countries, these may significantly deviate from the overall trend).

The proposed modeling framework uses a combination of the exponential growth model for the initial dynamics of the epidemic and a generic growth curve [8] to capture the observed patterns in the number of detected positive individuals. This phenomenological model is flexible, includes many special cases and thus allows selecting the effective parsimonious model fitting the observed data based on likelihood ratio tests [27,39] or information criteria such as the Akaike’s Information Criterion [40]. The effectiveness of this approach to phenomenological modeling has been demonstrated on COVID-19 data [5]. Our application on COVID-19 data from West Africa nevertheless showed that the logistic regression of recoveries and deaths in the identified positive individuals against time can lack fit, as measured by an asymptotic ${\chi }^2$ test on the residual deviance statistic. Nevertheless, these fits can be improved by adding explanatories (different from time, but related to available health facilities) in the logistic regression models. The deterministic SIQR model [22] considered for mechanistic modeling explicitly acknowledges the isolation of the detected positive individuals. It does not, however, include an exposed (E) state as in the SEIQR model [41]. The use of the SEIQR model may provide better insights on the effectiveness of control measures since most of the measures first impact the exposition of susceptible individuals. In general, the proposed modeling approach can be extended by considering more complex models such as the SEIQR and the SIDARTHE model [42] instead of the SIQR model considered herein.

Among interest quantities provided by the hybrid modeling framework, we have the epidemic latency period $t_o$ (the time from the appearance of the first infectious case in the population to the outbreak). For the West African region, the result indicates that the first imported COVID-19 case(s) in West Africa likely entered the region around 28 January–7 February 2020. To the best of our knowledge, this is the first estimate of this duration in the region. This epidemic latency period is much lower than the 40 days estimated for Italy [14]. This is in line with the relatively late arrival of the virus in the region, compared to the Asian and European continents, and the prevention and detection measures anticipated by many West African governments [31]. We obtained a basic reproduction number ($CI\left({\mathcal{R}}_o\right)=[2.60,2.69]$) higher than the estimate ($CI\left({\mathcal{R}}_o\right)=[1.84,1.87]$) obtained by [15]. Our estimate is, however, closer to country-specific estimates obtained for Nigeria ($CI\left({\mathcal{R}}_o\right)=[2.37,2.47]$) [43] and Ghana ($CI\left({\mathcal{R}}_o\right)=[1.99,3.37]$) [44].

During the early phase of the epidemic after the outbreak in West Africa, the detection and isolation of a fraction of infected individuals reduced the reproduction number from ${\mathcal{R}}_o$ to a control reproduction number of ${\widehat{\mathcal{R}}}_0=2.52$, i.e., about 5.26% decrease. We estimated the duration of this phase characterized by an exponential growth to be about one month after the outbreak. This implies that the control measures implemented by West African governments to limit the transmission of the disease were not effective on average before April 2020. Indeed, apart from measures taken to limit the importation of new positive individuals (travel bans), many actions to limit the local propagation of the disease were first implemented in late March 2020 [31] (e.g., curfew set up on 21 March in Burkina-Faso, on 23 March in Côte d’Ivoire, Mauritius and Senegal and on 26 March in Mali; city lockdown on 22 March in Ghana and on 29 March in Nigeria; isolation of the capital from the rest of the country in Côte d’Ivoire on 25 March 2020; and cordon sanitaire set up to isolate the south from the rest of the country on 30 March 2020 in Benin). Our results indicate that these measures started to impact the dynamics of the epidemic from early April 2020. However, the measures may have affected the transmission dynamics earlier, since the measures mainly limited the exposition of susceptible individuals to the disease.

After the exponential growth phase, the sub-exponential growth pattern allowed the epidemic to peak. The estimated peak time for the detected positive cases was around 15 July 2020, and close to the observed peak time (24 July 2020). This estimated date has a delay of about eight days with respect to the estimated peak time of new infections ($CI\left(t_{new}\right)=[126.18,136.11]$ days). This estimate is higher than the estimate ($CI\left(t_{new}\right)=\left[\mathrm{108,112}\right]$ days) obtained by [30]. These contrasting results may be related to the more realistic SIQR model considered in this work as compared to the simpler SIR model used by Honfo et al. [30] who ignored the quarantine-adjustment of the disease incidence [22]. On the contrary, the estimated maximum number of new infections ($CI\left({\dot{T}}_{max}\right)=[\mathrm{20,284.04},\mathrm{24,464.98}]$) agrees with the estimate ($CI\left({\dot{T}}_{max}\right)=[\mathrm{24,239},\mathrm{26,294}]$ new infections) obtained by Honfo et al. [30].

Our results show that the time-varying effective reproduction number has decayed over April–August 2020, reaching 1 on about 15 July 2020 and 0.66 at the end of the considered period (31 August 2020). Based on the modeled dynamics, the effective reproduction number likely reached its minimum value 0.61 around 29 September 2020. However, the reproduction number likely increased again to approach ${\mathcal{R}}_{\infty }=1$ in the long run. Overall, the various measures decided and enforced by different West African governments, against the first COVID-19 epidemic wave in the region, were able to contain the propagation of the disease (importation of new cases and local transmission) in five months.

However, the COVID-19 pandemic will remain an important issue for a long time, and local region’s endemic to the pathogen will likely appear in the long run. This is so because of the following factors: the re-opening of borders and airports in the region to limit the related economic feedback [45,46]; the relaxation of measures such as the ban of sport, political, cultural and religious gatherings [31,47]; and the natural evolution of the SARS-Cov-2 virus [48,49,50,51]. The limited resources and capacity of Sub-Saharan Africa countries in general [52,53,54] to immunize their population through vaccination will compound this threat in the region.

5. Conclusions

There are two common approaches to epidemiological modeling: phenomenological models and mechanistic models. This study proposes a hybrid framework which combines the two approaches, starting from fitting curves to observed data (confirmed positive cases, recoveries and deaths) and then providing an overall view of the epidemic dynamics by integrating the fitted curves into a compartmental model. The proposed approach allows estimating the delay between the appearance of the first infectious case in the population and the outbreak (“epidemic latency period”); the duration of period during which the epidemic growths exponentially; the basic and control reproduction numbers; and the peaks (time and size) in positive cases, active cases and new infections. An application to COVID-19 data from West Africa indicates that the hybrid modeling framework can be used to match effective control measures dictated by health policies with changes in the transmission dynamics of the studied disease.