key: cord-0802409-2894jei8
authors: Tovissodé, Chénangnon Frédéric; Lokonon, Bruno Enagnon; Glèlè Kakaï, Romain
title: On the use of growth models to understand epidemic outbreaks with application to COVID-19 data
date: 2020-10-20
journal: PLoS One
DOI: 10.1371/journal.pone.0240578
sha: bd89a1781a8678e0c8ea9a0a65eac87865b64213
doc_id: 802409
cord_uid: 2894jei8

The initial phase dynamics of an epidemic without containment measures is commonly well modelled using exponential growth models. However, in the presence of containment measures, the exponential model becomes less appropriate. Under the implementation of an isolation measure for detected infectives, we propose to model epidemic dynamics by fitting a flexible growth model curve to reported positive cases, and to infer the overall epidemic dynamics by introducing information on the detection/testing effort and recovery and death rates. The resulting modelling approach is close to the Susceptible-Infectious-Quarantined-Recovered model framework. We focused on predicting the peaks (time and size) in positive cases, active cases and new infections. We applied the approach to data from the COVID-19 outbreak in Italy. Fits on limited data before the observed peaks illustrate the ability of the flexible growth model to approach the estimates from the whole data.

COVID-19 is a pandemic caused by the new coronavirus strain SARS-nCOV2 which emerged from Wuhan, China [1, 2] . A total of 21 026 758 COVID-19 cases and 755 786 related deaths were reported across the world as at August 15, 2020 [3] . The worldwide social, as well as economic ravages by COVID-19 has immediately motivated the use of mathematical models to understand the course of the epidemic and plan for effective control strategies. These include, for instance, the SIR (Susceptible, Infectious, Recovered), SEIR (Susceptible, Exposed, Infectious, Recovered) and its variants, SIDR (Susceptible, Infectious, Recovered, Dead) and SIQR (Susceptible, Infectious, Quarantined, Recovered) models [4] [5] [6] [7] . These modelling approaches use mechanistic models which incorporate key physical laws or mechanisms involved in the dynamics of the population at risk and the pathogen [8] . A second class of approaches uses empirical phenomenological models which does not require specific knowledge on the physical laws or mechanisms that give rise to the observed epidemic data [9] , and was considered, for instance, by [10] and [11] to understand both short and long term dynamics of COVID-19. A new curve fitting-like approach, namely fractal interpolation [12, 13] was also proposed by [14] [15] [16] to account for the high noise and reporting bias in data from the COVID-19 pandemic. As generally is the case with dynamic biological systems [17, 18] , mathematical model development and adaptation are fundamental requirements to guide public health policies.

When facing an epidemic outbreak, public health officials are mostly interested in data driven, mathematically motivated, practical and computationally efficient approaches that can: i) generate estimates of key transmission parameters; ii) gain insight to the contribution of different transmission pathways; iii) assess the impact of control interventions (e.g. social distancing, test + isolation, vaccination campaigns); iv) optimize the impact of control strategies; and v) generate short and long-term forecasts [8] . In regard to the current COVID-19 outbreak, politics and public health officials are mostly worried about the ability of the disease to induce saturation of the health system, reducing the survival of patients, and even consulting for reasons different from the epidemic itself. High interest is thus currently given to accurate forecasting of the epidemic peak time and size, epidemic size and duration, as well as their sensitivity to control interventions in order to optimize the impact of control strategies.

An exponential-growth model is usually assumed to characterize the early phase of epidemics. But, this assumption can lead to failure to appropriately capture the profile of the epidemic growth, eventually giving rise to non-realistic epidemic forecasts [10, 19] . In an ultimate view to guide control interventions aiming to limit the spread of epidemics, with focus on the COVID-19 pandemic, this work considered a flexible growth curve fitting approach to understand the dynamics of epidemics. We used the generic growth model of [20] to model the course of reported positive cases and a binomial regression to model removals (recoveries and deaths). Thereafter, we inferred the overall dynamics of the epidemic, in terms of observables (reported cases, active/quarantined cases) and unobservables (new infections, lost cases), and predicted interest quantities such as the peak (time and size) in reported cases, active cases and new infections. The performance of the approach was assessed through an application to daily case reporting data from Italy, which has virtually completed a whole COVID-19 outbreak wave, thus offering the possibility to compare predicted outputs to real events.

We used a growth curve approach for modeling the course of an epidemic along time. We followed [8, 10, 21] who, among others, used growth models to forecast epidemic dynamics.

Let C t denote the size of the detected infected population at time t, i.e. the cumulative number of infected, identified and isolated individuals. We assumed for convenience that C t is continuous and denote _ C t its first derivative with respect to t. Also let I t be the true size of infectives at t, related to C t through

where δ t 2 (0, 1] is the detection rate which is closely related to the testing effort (number of tests, tracing of contact persons of identified cases and targeting exposed people) and is assumed at least twice differentiable with respect to t. We ressorted to the generic growth model of [20] for the identified positive cases:

with u t = [1 + νωρ(t − τ)] −1/ρ . In Eq (2), K > 0 is the ultimate epidemic size (detected), ω > 0 is the "intrinsic" growth constant, ν and ρ are powers (ν > 0 and −1 < ρ < ν −1 ) characterizing respectively the rates of change with respect to the initial size C 0 = δ 0 I 0 (number of cases detected at time t = 0) and the ultimate size K, and τ is a constant of integration, determined by the initial conditions of the epidemic and implicitly the detection rate δ 0 through C 0 = K[1 + (1 − νωρτ) −1/ρ ] −1/ν for ρ 6 ¼ 0 and C 0 = K(1 + e νωτ ) −1/ν for ρ = 0. The growth model in Eq (2) is quite flexible to handle various shapes of epidemic dynamics. Indeed, if K ! 1 and νρ ! 0, Eq (2) specializes to the exponential growth model

where ω is the exponential growth rate. Apart from Eq (3), other special or limiting cases of Eq (2) include the hyper-Gompertz (ν ! 0 while ων 1+ρ is constant) and the Gompertz (ν ! 0, ρ ! 0 while ων is constant), the Bertalanffy-Richards (ρ ! 0), the hyper-logistic (ν = 1) and the logistic (ν = 1 and ρ ! 0) growth models [20] . From Eq (2), the observed epidemic incidence _ C t is given by

In order to ensure the restriction −1 < ρ < ν −1 , we set r ¼ r 0 nþ1 n À 1 with ρ 0 2 (0, 1) free of ν.

The number A t of detected and active cases along an epidemic outbreak is of high interest for public health officials. Indeed, A t must be kept under the carrying capacity of the health system to avoid overload and disrupture. The derivative _ A t of the detected and active cases satisfies

where R t = α t A t denotes the number of removed and permanently immune (mortality and recovery) at time t, and α t is the unit time removal probability, i.e. the odds to have an outcome (recovery or death), averaged over the active cases. Eq (5) fits in the SIQR (Susceptible, Infectious, Quarantined, Recovered) model framework [22] with the detected active cases referred to as "quarantined" and the strong assumption that α t is constant along the epidemic outbreak (see the third equation in system (6) in [22] ). The removal probability can more generally be given the logistic form a t ¼ e Z t 1þe Z t with Z t ¼ X > t β þ kt where X t = (X t1 , X t2 , � � �, X tq ) > is a vector of q covariates (known constants) and β is the q vector of associated effects, and κ determines the change in the log-odds ratio for having an outcome per unit time. These changes in α t can be due to an improvement in the health care system during the epidemic outbreak (increase in recovery ratio) or a deterioration of the health care system for infected individuals (increase in mortality ratio due to the outbreak). The general solution of the differential Eq (5) turns to have the form

where A 0 is the number of active cases at time t = 0. Indeed, when κ = 0, taking the first deriva-

C s e a s s ds�ðÀ a t Þe À a t t resulting in _ A t ¼ _ C t À a t A t which is the Eq (5). For t = 0, the integral in Eq (6) vanishes, resulting as expected in A t = A 0 since e À a t t¼1 . When κ 6 ¼ 0, the first derivative of A t is _

Eq (5). Here, for t = 0, e Z t ¼ e X > t β so that A t = A 0 . There are no general closed form solutions for the integrals in Eq (6), unless _ C t and α t are purposely chosen as functions of time to simplify the integral. A t can, however, be obtained in practice from Eq (6) using a numerical integration routine such as the function integrate in R freeware [23] or the function integral of Matlab [24] . Nevertheless, to circumvent this issue during estimation under the generic growth model in Eq (2), we discretized the active cases A t by assuming a binomial removal process R t conditional on the detected unit time new cases Y t as

where BIN(n, α) denotes a binomial distribution with n trials and success probability α and Y t is a non-negative process with expectation l t ¼ _ C t . Clearly, the bivariate process {A t , R t } defined by Eqs (8) and (7) is not stationary. However, since Y t � 0 and _ C t ! 0 as t ! 1, we have Y t ! 0 in distribution as t ! 1, and if the removal probability α t does not approach zero as t ! 1, then A t ! 0 as t ! 1.

The epidemic peak is an important event in the disease dynamic and can be estimated for a better management of the epidemic. An epidemic described by the exponential growth model in Eq (3), (K ! 1 or νρ ! 0) does not peak. Otherwise, the peak in the detected number of infected individuals corresponds to the maximum of the incidence rate _ C t . This maximum is

Solving € C t ¼ 0 for t using Eq (9) yields the peak time

on replacing r ¼ r 0 nþ1 n À 1. Inserting t p in Eq (1) and denoting u p ¼ n 1þr 1À rn gives the peak

At the peak in detected cases, the cumulative number of detected cases is

An important interest in modelling the epidemic incidence is the derivation of quantities related to the overall dynamics of the epidemic, in both detected and undetected cases. Total cases: Detected and losts. Let us denote S t the cumulative number of cases from the epidemic outbreak to t, and let _ S t be the first derivative of S t . We also introduce Λ t , the cumulative number of lost cases (with first derivative _ L t ), i.e. people who were infected, undetected, and removed from infectives (mortality and recovery).

The size of the lost cases is determined by the unit time removal rate π t 2 (0, 1) from undetected infectives (π t is an average over all infectives, i.e. irrespective of the time since infection onset). The lost rate π t which is assumed at least twice differentiable with respect to t, depends on various factors like the disease related mortality, the average infection duration, the natural proportion of asymptomatics within infectives, and the existence and the use of medicines that may reduce symptoms (induced asymptomatics). It is worthwhile noticing that π t can be estimated from the removal rate α t in the detected cases, taking into account various factors that may induce difference between the two rates. For instance, since the undetected cases include asymptomatics, disease related mortality may be lower and recovery rate higher in undetected as compared to detected cases. However, efficiency of the health care system in treating identified and isolated cases can reduce mortality thereby reducing α t , but also improve recovery thereby increasing α t .

With the above notations, the lost cases count Λ t satisfies the differential equation,

whereas the cumulative number of cases S t is given on

The factor υ t represents at time t the proportion of infectives who will potentially continue to spread the epidemic after adequate contacts (i.e. contacts sufficient for transmission) with susceptibles. In other words, the number of undetected currently infectives is (1), the infectives I t and its first derivative with respect to time _ I t are given for t � 0 by

where _ d t is the first derivative of the detection rate δ t with respect to t. Straightforward algebraic operations then give the number of new cases and the cumulative number of cases as

p t the first derivative of the lost rate π t , and the cumulative number of lost cases Λ t is given for t � 0 by

with S 0 the cumulative number of all cases until the first detection date t = 0. The total size of the epidemic is

Under the Turner's growth model,

Let us assume a constant detection rate δ t = δ closely related to detection effort but also to the average duration from infection to recovery or death of non-isolated cases. Assuming in addition a constant lost rate (π t = π), we have

, and the new cases _ S t and its accumulation S t , as well as the lost cases Λ t simplify to

The total epidemic size is here

Epidemic peak. At the time t p of the peak of reported cases ( € C t ¼ 0) under constant detection and lost rates, the new infectives is _ S p ¼ ½1 þ pðd À 1 À 1Þ� _ C p with _ C p given in Eq (11) . This, however, corresponds to the peak in the overall new cases _ S t only under the unrealistic assumption δ = 1. The peak of new infections occurs when the second derivative € S t of S t with respect to t vanishes ( € S t ¼ 0). We have from Eq (16)

where C ⃛ t (the third derivative of C t with respect to t) and C t are given by

with z t ¼ u t 1þu t , and € p t and € d t the second derivatives of respectively π t and δ t with respect to t.

The peak time and value depend on the particular forms of δ t and π t as functions of time.

Here, we restrict the attention to the simple situation with constant positive detection and lost rates (δ t = δ with δ 2 (0, 1) and π t = π with π 2 (0, 1)) where

It appears that the peak of new infections occurs before the time t p of the peak in detected cases. Indeed, at t = t p , we have

_ S t is already in its descending phase. The expression Eq (25) indicates that at the time t P of the peak of new infections, € C t is equal to € C P ¼ À zC ⃛ P where z ¼ ð1À pÞð1À dÞ pþdð1À pÞ and C ⃛ P is given by Eq (23) with t = t P . The lower z, the lower j € C P j, and the lower the difference t p − t P (delay of the observed peak). Differentiating z with respect to δ gives @z @d ¼ À 1À p ½pþdð1À pÞ� 2 < 0, hence the higher δ, the lower the delay between the observed peak time and the time of the peak in new infections. Using Eqs (9) and (23), € S t becomes

which does not have a closed form root. The root t P can, however, be obtained using root finding numerical routines such as the R function uniroot or the Matlab function fzero. Afterwards, the peak _ S p size (the maximum number of new infections) is obtained using Eq (19) .

Let us consider a record of new confirmed infected cases Y 1 , Y 2 , � � �, Y n , active cases A 0 , A 1 , � � �, A n−1 , removed cases R 1 , R 2 , � � �, R n (available from Eq (8) as R t = Y t − A t + A t−1 ) and the associated vectors of covariates X 1 , X 2 , � � �, X n at n time points. The parameters K, ω, ν, ρ 0 , τ, and κ can be estimated using maximum likelihood (ML) by assigning to each Y t an appropriate statistical distribution with expectation l t ¼ _ C t and a dispersion parameter σ > 0, and probability density function (pdf) or probability mass function (pmf) f(Y t |θ) where θ = (K, ω, ν, ρ 0 , τ, κ, β > , σ) > . We subsequently considered inference under log-normal and negative binomial distributions.

Log-normal model. Epidemic incidence case data are generally fitted through non-linear least squares applied at logarithmic scale [19, 25, 26] . To deal with zero incidence cases, the logarithmic transform is usually applied on the shifted cases Y t + 1. Mimicking this procedure in a likelihood inference framework, we consider a log-normal distribution assumption for the shifted incidence cases, i.e. Y t + 1 * LN(λ t + 1, σ). The pdf of Y t , adapted from [27] , reads

ffi ffi ffi ffi ffiffi 2p p exp À 1 2

so that Y t has expectation E[Y t ] = λ t and variance Var½Y t � ¼ ðl t þ 1Þ 2 ðe s 2 À 1Þ.

Negative binomial model. Since incidence cases are counts, Y t can be assumed to follow the negative binomial distribution, i.e. Y t * NB(λ t , σ) with pmf

The incidence case Y t then has expectation E[Y t ] = λ t and variance Var[Y t ] = λ t (1 + σλ t ). Likelihood inference. Based on the information {Y t , R t } for t = 1, 2, � � �, n, the conditional log-likelihood of the parameter θ given A 0 is

Þa R t t ð1 À a t Þ A tÀ 1 þN t À R t is the binomial probability mass function for R t .

The function ℓ(�) can be maximized to obtain the maximum likelihood estimateθ of θ using an optimization routine such as the function optim in R or the function fminsearch of Matlab.

Let H(θ) the hessian matrix of ℓ(θ) and define the covariance S(θ) = −[H(θ)] −1 . The large sample distribution (i.e. for n ! 1) of the maximum likelihood estimator is multivariate normal with meanθ and covariance matrixŜ ¼ SðθÞ.

The data. In order to test the reliability of the Turner's growth model in predicting the dynamics of an epidemic, we used data from one of the countries which had completed a whole COVID-19 outbreak wave. The daily case reporting data in Italy was obtained from https://github.com/CSSEGISandData/COVID-19/tree/master/csse_covid_19_data/csse_ covid_19_time_series. We used only the confirmed data (2020-02-20 to 2020-07-11) accessed on 2020-07-28, discarding the latest data subject to possible reporting delay, as indicated by the Istituto Superiore di Sanità (ISS) at https://www.epicentro.iss.it/en/coronavirus/sars-cov-2-dashboard.

Data analysis. All analyses were performed in R [23] . We fitted the Turner's growth model curve to the whole Italian data. Both the log-normal and the negative binomial distributions were used, and the fit with the lowest root mean square error (RMSE) computed for the daily new positive cases was selected as the best. Then, we derived peak statistics (time and size) for daily new reported cases and active cases. We also inferred the daily new infections from assuming constant detection and lost rates and estimated its peak (time and size). The detection (δ = 0.033/day) and the lost rates (π = 0.1/day) for Italy were obtained from [7] . These rates follow from assumptions that the average time of duration from infection to recovery or death of non-isolated cases is 10 days (hence π = 0.1/day) and that during this detection window, 1/3 of infectives are tested positives (hence δ = 0.033/day).

To assess the ability of the model in predicting the peak of the new positive cases in countries which have not yet reached the peak, we retrospectively fitted the model to the Italian data before the observed peak (day 29 after the notification of first case), using data of the first two weeks, and then data of the first three weeks. For these analyses with limited data, we fitted the full Turner's growth model to the positive cases, but also its special cases, namely the hyper-Gompertz (ν ! 0 while ων 1+ρ is constant), the Gompertz (ν ! 0, ρ ! 0 while ων is constant), the Bertalanffy-Richards (ρ ! 0), the hyper-logistic (ν = 1) and the logistic (ν = 1 and ρ ! 0) models using the log-normal distribution for the daily counts. We then computed the Akaike's Information Criterion (AIC) defined as AIC ¼ À 2' þ 2N p with' the maximized loglikelihood and N p the number of parameters in a fitted model. Finally, we retained and presented the best fit (lowest AIC value). Table 1 shows parameter estimates using the whole Italian COVID-19 daily case reporting data from 2020-02-20 to 2020-07-11, with standard errors and 95% confidence intervals. The log-normal distribution based fit recorded the lowest RMSE, and was thus retained for subsequent analyses. The confidence bounds for the parameter ρ Notes: RMSE = root mean square error; β and κ define the daily removal rate from detected cases as a t ¼ e bþkt 1þe bþkt ; σ is the log-normal/negative binomial distribution scale parameter (see the pdf in Eq (27) and pmf in Eq (28)).

https://doi.org/10.1371/journal.pone.0240578.t001 ), it appears that the daily removal rate (recoveries and deaths) averagedâ 0 ¼ 1:8% in the very early phase of the epidemic (t � 0 day). Then, from the estimate of κ (k ¼ 0:0076 with CI(κ) = [0.0075, 0.0076]), it appears that the removal rate increased with time, i.e. the probability for an active case to recover or die within a day increased on average by 5.5% over a week . Fig 1(C) displays the active cases and the corresponding fitted curve using the removal probability along with the fitted Eq (30). The active cases were predicted to peak on day 56 (t a ¼ 55: The daily new infections inferred from assuming a constant detection rate (δ = 0.033/day) and a constant lost rate (π = 0.1/day) is depicted on Fig 1(D) 

The AICs of the retrospective fits of Tuners's growth model and its special cases to the Italian COVID-19 data of the first two weeks and the first three weeks are presented in Table 3 . It can be observed that the best fits correspond to the hyper-logistic growth model for both data of the first two weeks (AIC = 483.03) and data of the first three weeks (AIC = 863.58). Although parsimony indicated the hyper-logistic model fits as the best, the differences Δ AIC in AIC with respect to the full Turner's growth model fit were mild (|Δ AIC | < 2). Table 4 shows the estimate of the hyper-logistic growth model parameters for the two shorted datasets. It appears that the estimates of the intrinsic growth parameter ω increased slightly with data availability fromô ¼ 0:05 (CI ω = [0.04, 0.07]) using the data of the first two weeks, toô ¼ 0:07 (CI ω = [0.06, 0.08]) using the data of the first three weeks and toô ¼ 0:09 (CI ω = [0.07, 0.11]) using the whole dataset from Italy. The estimates of the peak time and size from the two shorted datasets are shown in Table 4 . The forecast of the peak time from the data of the first two weeks was day 44 (t p ¼ 43 (Table 5 ).

This work proposes the use of a flexible growth model to model case reporting data from an epidemic outbreak with containment measures including at least isolation of individuals tested positive. The generic growth model of [20] offers a flexible framework with the possibility to recover many special growth models such as the common exponential and the logistic growth models, the hyper-logistic, the hyper-Gompertz, the Gompertz and the Bertalanffy-Richards growth models. Since the special models are all nested within the generic model framework, the most appropriate model can be identified using information criteria such as the Akaike's Information Criterion (AIC), but a likelihood ratio test [28] can also be conducted for models with different number of free parameters. Where additional information can be obtained on the ability to detect infective individuals, the proposed framework allows to include this information so as to infer on the dynamics of the epidemic beyond the identified (positive) cases, without ressorting to mechanistic/compartmental models. Nevertheless, we considered a constant (average) detection rate whereas the detection rate obviously changes over the epidemic course in terms of the detection effort (number of tests, tracing of contact persons). From our application to the COVID-19 outbreak data in Italy, the hyper-logistic model is the most appropriate model for the dataset. It appears that the modelling approach can predict the dynamics of an epidemic using data from first few days of an outbreak, at least in this example. Indeed, the predicted peak time (and size) for the positive cases (using only the first two/three weeks data) overestimates (and underestimates) the observed peak time (and size). However, the biases can be attributed, for instance, to the increase in the testing effort and isolation (and the subsequent decrease in the growth rate) in Italy where only about 3 762 tests/day were performed in the first three weeks from 2020-02-20, and about 21 248 tests/day were performed in the subsequent three weeks. Our estimate of the ratio of the number of infectives to the number of active cases averaged 22.95 in the first week of the outbreak, within the range [5, 25] obtained by [7] using the SIQR model. Our proposal thus offers a valid alternative to mechanistic models, for instance, the piecewise exponential growth used by [7] within the SIQR model framework on the Italian early outbreak data.

In a very limited data situation, we suggest a further reduction of the number of model parameters to be estimated. Indeed, since the parameter τ in the growth model in Eq (2) is a constant of integration determined by the initial conditions of the epidemic, it can be expressed in terms of other parameters and the number of cases C 0 detected at time t = 0 as

h i À r n o for ρ 6 ¼ 0 and τ = log((K/C 0 ) ν − 1)/(νω) for ρ = 0. Consideration of a procedure where τ is not estimated as a free parameter may lead to parsimony, with inference conditional on the number of individuals tested positive at time t = 0. Inference on the effective reproduction number and the sensitivity of the epidemic dynamics to the containment measures under the generic growth model framework is considered for future work. 

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

The COVID-19 epidemic. Tropical medicine & international health

WHO. Coronavirus disease 2019 (COVID-19): situation report

Data-based analysis, modelling and forecasting of the COVID-19 outbreak

Can the COVID-19 epidemic be controlled on the basis of daily test reports?

Early dynamics of transmission and control of COVID-19: a mathematical modelling study. The lancet infectious diseases

Quantifying undetected COVID-19 cases and effects of containment measures in Italy

Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecasts. Infectious Disease Modelling

Mathematical models to characterize early epidemic growth: A review

Modeling the Covid-19 Epidemic using

A Poisson autoregressive model to understand COVID-19 contagion dynamics. SSRN ePrint

Primality, fractality, and image analysis

Entropy and fractal antennas

An analysis of COVID-19 spread based on fractal interpolation and fractal dimension

Some fractal thoughts about the COVID-19 infection outbreak

A fractal kinetics SI model can explain the dynamics of COVID-19 epidemics

Immune system network and cancer vaccine

Persistence analysis in a Kolmogorov-type model for cancer-immune system competition

Is it growing exponentially fast?-impact of assuming exponential growth for characterizing and forecasting epidemics with initial near-exponential growth dynamics. Infectious disease modelling

A theory of growth

Real-time forecasting of epidemic trajectories using computational dynamic ensembles

Effects of quarantine in six endemic models for infectious diseases

R: A Language and Environment for Statistical Computing

Comparative estimation of the reproduction number for pandemic influenza from daily case notification data

A generalized-growth model to characterize the early ascending phase of infectious disease outbreaks

Log-normal Distributions across the

The large-sample distribution of the likelihood ratio for testing composite hypotheses. The annals of mathematical statistics

The authors are grateful to an anonymous reviewer for drawing their attention on models based on fractal-wavelet for modelling the COVID-19 pandemic. RGK acknowledges the support from the African German Network of Excellence in Sciences (AGNES).

Conceptualization: Chénangnon Frédéric Tovissodé, Romain Glèlè Kakaï.