“All models are wrong but some are useful” / “Todos os modelos estão errados, mas alguns são úteis” (George E. P. Box)

“Prediction is very difficult, especially if it’s about the future” / “Previsão é muito difícil, especialmente sobre o futuro” (Niels Bohr)

As citações anteriores, atribuídas ao estatístico Box e ao físico Bohr, retratam a grande dificuldade em modelar e prever a maioria dos fenômenos naturais, biológicos e sociais. No contexto particular da pandemia do COVID-19 muitos modelos têm sido resgatados e novos modelos têm sido propostos para tentar prever o avanço da doença. Ainda assim, não existe nenhum modelo à prova de falhas para descrever a evolução dos casos de COVID-19.

Nesse contexto, apresentamos uma descrição um pouco mais detalhada de algumas possibilidades para descrever o comportamento de doenças infecciosas usando modelos matemáticos.

Modelo exponencial

A grande maioria das epidemias apresenta um crescimento exponencial em seu início [1]. Essa característica reflete o fato de que, inicialmente, nenhum indivíduo está imune à doença e aqueles que contraem o vírus geralmente o transmitem para mais de uma pessoa.

A figura abaixo mostra o número acumulado de pessoas infectadas em função do tempo (círculos em preto) para os 20 primeiros dias após o primeiro caso em Maringá. Nesse gráfico, a linha tracejada mostra o modelo exponencial ajustado aos dados. Matematicamente, esse modelo é expresso por

$$ I(t) = I_0~\exp({r~t}),$$

na qual $I$ representa o número de infectados, $t$ o tempo em dias, $I_0$ é uma constante que representa o número inicial de infectados [$I(t=0)=I_0$] e $r$ representa a taxa de crescimento do número de casos. A taxa $r$ indica o quão rápido o número de casos está crescendo: quanto maior o valor de $r$, mais rapidamente o número de casos aumenta com o passar do tempo.

Além disso, a área sombreada em azul mostra a região de incerteza do modelo. Podemos perceber que a descrição exponencial ajusta-se bem ao início da curva do número de casos. Entretanto, esse modelo passou a superestimar o número de infectados após esse período inicial. Isso pode ser um indicativo de desaceleração da epidemia, mas também pode ser um artefato causado pela escassez de testes da doença.

Definição do tempo de dobra

Uma outra maneira de escrever o modelo exponencial para o número de casos é

$$ I(t) = I_0~2^{~t/t_2}, $$

sendo o parâmetro $t_2$ conhecido como tempo de dobra ou período de dobra. Como sugere o nome, o tempo de dobra $t_2$ indica o tempo necessário para que o número de casos dobre. Assim, se o tempo vai de $t$ para $(t+t_2)$, o número de casos aumenta de $I(t)$ para $2\times I(t)$. De modo similar à taxa de crescimento $r$, o tempo de dobra também indica quão rápido o número de casos está crescendo. Nesse caso, quanto menor o valor do tempo de dobra, mais rapidamente o número de casos aumenta no tempo. Vale notar que o tempo de dobra $t_2$ e a taxa de crescimento $r$ estão relacionados via $r = \ln 2 / t_2$. A terceira figura da nossa página sobre a evolução da COVID-19 em Maringá mostra a evolução da estimativa para o tempo de dobra em nossa cidade.

Modelo logístico

Uma outra possibilidade para descrever o número de infectados em função do tempo é usar o modelo logístico. Uma versão alterada desse modelo, conhecido como modelo generalizado de Richards, foi recentemente utilizado para modelar o número de casos de COVID-19 em províncias da China, em países como Japão, Coreia do Sul e Irã e na Europa como um todo [2]. No caso mais simples, temos que o modelo logístico é descrito matematicamente por uma equação diferencial

$$ \frac{dI}{dt} = r I (1-I/K) $$ na qual $r$ representa a taxa de crescimento e $K$ é uma outra constante denominada capacidade de carga. A solução da equação anterior é normalmente denominada função logística e pode ser escrita como

$$ I(t) = \frac{K}{1+\frac{K-I_0}{I_0} \exp(rt)}, $$ na qual $I_0$ representa o número inicial de infectados.

O modelo logístico pode ser entendido como uma generalização do modelo exponencial, que é recuperado quando o valor de $K$ fica muito grande ($K\to\infty$). Para valores finitos de $K$, o modelo logístico apresenta um crescimento exponencial inicial (tal como o que observamos na figura anterior) que passa a saturar após algum tempo. Após um longo tempo, esse modelo prevê que o valor de $I$ se aproxima da constante $K$, ou seja, $I(t\to\infty)\to K$. Desse modo, percebemos que o modelo logístico é uma descrição um pouco mais realista para a evolução da epidemia, uma vez que o número de casos não pode aumentar para sempre numa população com um número finito de pessoas.

A figura abaixo mostra o ajuste do modelo logístico para o número de casos de COVID-19 em Maringá considerando os primeiros 30 dias após o primeiro caso na cidade. Nesse gráfico, os círculos em preto indicam os dados reais e a linha tracejada mostra o comportamento do modelo logístico. Além disso, a área sombreada em azul mostra a região de incerteza do modelo. Notamos que o modelo logístico descreve o número de casos razoavelmente bem desde o início da epidemia até o dia 18/04/2020. Entretanto, após essa data o modelo logístico passou a subestimar o número casos. Por volta de 21/04/2020 o número de casos apresentou um segundo ponto de inflexão, o qual não pode ser descrito pelo modelo logístico. Esse segundo ponto de inflexão indica uma possível retomada no crescimento de casos de COVID-19 na cidade.

Modelo bi-logístico

Como já mencionamos, a curva do número de casos acumulados ao longo do tempo de COVID-19 em Maringá apresenta mais de um ponto de inflexão. Essa característica faz com que o modelo logístico usual (descrito na seção anterior) seja incapaz de descrever o número de casos de modo adequado. Devido a esse comportamento e inspirados em um trabalho de modelagem do espalhamento da COVID-19 na Itália[3], resolvemos usar uma generalização do modelo logístico conhecida como modelo multi-logístico (mais particularmente o modelo bi-logístico). De modo qualitativo, modelos multi-logísticos descrevem processos de crescimento com mais de um regime de crescimento e saturação. No contexto particular da pandemia de COVID-19, esses diversos regimes podem estar associados à adoção de medidas de distanciamento social (ou mesmo lockdown) seguidos de suas flexibilizações ou também ao surgimento de novos focos de propagação da doença à medida que o tempo passa.

Do ponto de vista matemático, modelos multi-logísticos podem ser expressos pela soma de funções logísticas defasadas no tempo, nas quais os tempos de defasagens representam a separação temporal entre cada um dos regimes de crescimento. No modelo bi-logístico, existem dois regimes de crescimento que podem ser escritos pela soma $$ \begin{align} I(t) &= \frac{K_0}{1+\frac{K_0-I_0}{I_0} \exp(r_0 t)} \\
&+ \frac{K_1}{1+\frac{K_1-I_0}{I_1} \exp(r_1 (t-\tau_2))}. \end{align} $$ Note que essa expressão nada mais é que a soma de duas funções logísticas com seus respectivos parâmetros. Além disso, $\tau_2$ representa a defasagem temporal entre os dois regimes de crescimento logístico. Esse modelo representa uma descrição fenomenológica mais geral para processos de crescimento, uma vez que permite descrever uma curva de crescimento com mais de um ponto de inflexão, tal qual observamos nos dados de Maringá.

Verificamos que o modelo bi-logístico é uma descrição razoável para a evolução do número de casos acumulados de COVID-19 em Maringá ao considerarmos os primeiros 55 dias dias após o primeiro caso. Porém, após essa data, o modelo passou a não representar bem a evolução do número casos da doença. Como mencionado anteriormente, existem dificuldades para se encontrar um modelo que se ajuste à evolução da epidemia. O espalhamento da doença ocorre num ambiente não controlado e as diversas ações que vêm sendo tomadas pelos governantes têm um impacto não-trivial no espalhamento da COVID-19. Normalmente, é muito difícil levar esses diversos fatores em conta nos modelos. Também, sabemos que existe uma grande subnotificação dos casos da COVID-19 no mundo, no Brasil e, certamente, em Maringá, o que também afeta qualquer tentativa de prever o comportamento futuro da pandemia. A figura abaixo mostra o ajuste do modelo bi-logístico para o número de casos de COVID-19 em Maringá do início da epidemia até o dia 11/05/2020. Nesse gráfico, os círculos em preto indicam os dados reais e a linha tracejada mostra o comportamento descrito pelo modelo bi-logístico. Além disso, a área sombreada em azul mostra a região de incerteza do modelo.

Modelo epidêmico SIR

O modelo epidêmico SIR é um dos modelos epidêmicos mais simples, mas é capaz de descrever diversas características encontradas em epidemias. Esse modelo assume que a população está homogeneamente misturada (todos podem interagir com todos) e que cada indivíduo da população está Suscetível (ou saudável), Infectado ou Recuperado (imune) da doença (daí a sigla SIR). Sua dinâmica assume que indivíduos suscetíveis podem se tornar infectados ao entrar em contato com indivíduos infectados ($S\to I$ ), e que indivíduos infectados se recuperam ($I\to R$ ) após um período de tempo, tornando-se potencialmente imunes. Matematicamente, podemos escrever a dinâmica desse modelo por um conjunto de três equações diferenciais $$ \begin{align} \frac{d S(t)}{dt} &= -\beta S(t) I(t)/N,\\
\frac{d I(t)}{dt} &= \beta S(t) I(t)/N - \gamma I(t),\\
\frac{d R(t)}{dt} &= \gamma I(t), \end{align} $$ no qual $S(t)$ representa o número de indivíduos suscetíveis no tempo $t$, $I(t)$ é o número de infectados em $t$, $R(t)$ é o número de recuperados no tempo $t$ e $N = S(t)+I(t)+R(t)$ representa a população total (uma constante). O modelo possui dois parâmetros: $\beta$ é a taxa de contato e $1/\gamma$ é o período infeccioso. A taxa de contato $\beta$ é uma propriedade da população que está sendo infectada e representa o número médio de contatos sociais por unidade de tempo que podem transmitir a doença. Já o período infeccioso $1/\gamma$ é uma propriedade fisiológica da doença e representa o tempo médio de duração da infecção.

Número de reprodução básico e tamanho da epidemia

Ainda nesse mesmo contexto, é comum definir o chamado número de reprodução básico $\mathcal{R}$ como sendo a razão entre $\beta$ e $\gamma$, ou seja, $$ \mathcal{R} = \frac{\beta}{\gamma}. $$ O número de reprodução básico representa o número médio de novas infecções provocadas por um indivíduo infectado na população na fase inicial da epidemia. Para a COVID-19, o valor de $\mathcal{R}$ tem sido estimado entre $2.2$ e $6.5$ [4], um valor muito maior do que o usualmente reportado para gripe [5] (entre 0.9 e 2.1). O valor de $\mathcal{R}$ é muito importante para a dinâmica do modelo SIR, particularmente quando estamos na fase inicial da epidemia e a maioria dos indivíduos é suscetível à doença (matematicamente, $S/N\approx 1$). Nesse caso, podemos aproximar a evolução do número de infectados pela equação $$ \frac{dI}{dt} \approx \gamma(\mathcal{R}-1)I, $$ cuja solução é uma função exponencial $$ I(t) \approx I_0 \exp(\gamma(\mathcal{R}-1) t). $$ Dessa última expressão, notamos que o número de casos cresce exponencialmente se $\mathcal{R}>1$. Por outro lado, o número de casos decresce exponencialmente quando $\mathcal{R}<1$. O valor $\mathcal{R}=1$ é normalmente chamado de limiar epidêmico.

Usando o valor de $\mathcal{R}$ também podemos calcular o tamanho da epidemia, ou seja, o número de indivíduos suscetíveis após um longo período de tempo ($t\to\infty$). Essa quantidade pode ser estimada resolvendo a equação [6] $$ \ln \frac{S_\infty}{S_0} = -\mathcal{R} (1-S_\infty/N), $$ na qual $S_\infty$ representa o total de indivíduos suscetíveis no final da epidemia (as pessoa que não foram infectadas). Subtraindo $S_\infty$ do total de indivíduos na população ($N$) podemos estimar a fração esperada de pessoas que serão infectadas até o final da epidemia. Por fim, podemos calcular também o tempo necessário para que toda essa fração da população seja infectada. Nessa data, denominada data de saturação da epidemia, espera-se que o número de novas infecções seja muito pequeno. Duas figuras abaixo apresentam estimativas para a porcentagem da população que será infectada e para a data de saturação da epidemia.

Novamente, cabe ressaltar que essas estimativas sofrem de problemas ligados à provável subnotificação de casos da COVID-19, além da simplificação exagerada imposta pela dinâmica do modelo SIR que usamos.

Número de reprodução instantâneo

Uma outra opção para o cálculo do número de reprodução é a escolha de abordagens que levam em conta o perfil de infectividade dos contaminados, também conhecidos como modelos time-since-infection [9]. A partir deles, é possível calcular o número de reprodução instantaneamente no tempo $t$, que apresenta a mesma interpretação do número de reprodução básico, porém tem validade para qualquer momento da epidemia. Para isso, modelamos a transmissão de COVID-19 como um processo de Poisson por meio de uma abordagem Bayesiana. O número de reprodução instantâneo pode ser descrito como $$ \mathcal{R}(t) = \dfrac{I(t)}{\sum_{s=1}^{\tau} I(t-s)w(s)}\ , $$ em que $w(s)$ é o perfil de infectividade e $\tau$ é a janela temporal considerada para estimar o $\mathcal{R}(t)$, considerando que o número de reprodução instantâneo permanece constante no intervalo $[t-\tau+1, t]$. O perfil de infectividade $w(s)$ é uma distribuição de probabilidade que depende do tempo desde a infecção por COVID-19. Dessa forma, maiores valores de $w(s)$ indicam que os infectados $I(t-s)$ estão mais suscetíveis a transmitir a doença nesse período. A distribuição $w(s)$ depende apenas de fatores biológicos individuais como características do vírus e severidade dos sintomas [10]. Dessa maneira, o número de reprodução instantâneo $\mathcal{R}(t)$ como calculado na equação acima equivale ao número médio de casos gerados pelos infectados na janela de $\tau$ dias anteriores. A verossimilhança para $I(t)$ pode ser escrita como uma distribuição de Poisson, isto é, $$ \mathcal{L}[I(t)| I(t-\tau+1),\dots,I(t-1), w, \mathcal{R}(t)] = \dfrac{[\mathcal{R}(t)\Lambda(t)]^{I(t)}e^{\mathcal{R}(t)\Lambda(t)}}{I(t)!}\ , $$ em que $$ \Lambda(t)=\sum_{s=1}^{\tau}I(t-s)w(s)\ . $$ Dessa forma, supondo independência estatística entre os períodos, a verossimilhança dentro da janela $\tau$ pode ser escrita como $$ \mathcal{L}[I(t-\tau+1),\dots,I(t)| w, \mathcal{R}(t)] = \prod_{s=t-\tau+1}^{t} \dfrac{[\mathcal{R}(t)\Lambda(s)]^{I(s)}e^{-\mathcal{R}(t)\Lambda(s)}}{I(s)}\ . $$

Usando o teorema de Bayes, podemos calcular a distribuição a posteriori para estimar $\mathcal{R}(t)$, isto é,

$$ P[\mathcal{R}(t),I(t-\tau+1),\dots,I(t)| w ] \propto \mathcal{L}[I(t-\tau+1),\dots,I(t)| \mathcal{R}(t), w ] P[\mathcal{R(t)}]\ . $$

Seguindo o procedimento proposto pelos sites epiforecasts.io e covid19br, consideramos que o perfil de infectividade é uma distribuição Gamma com média de 4,8 dias e desvio padrão de 2,3 dias, valores estimados em [7]. Além disso, atribuímos incerteza ao perfil de infectividade amostrando a média e o desvio padrão de distribuições normais truncadas centradas nos valores descritos anteriormente para um total de 500 repetições. Consideramos que $\mathcal{R}(t)$ apresenta uma distribuição Gamma como priori, com média 2,6 e desvio padrão 2, o que leva em consideração estimativas do número de reprodução básico em Wuhan no começo da epidemia [11]. Por fim, para a janela temporal $\tau$, adotamos o tempo $\tau=7$ dias uma vez que a série temporal de casos por COVID-19 aparenta uma periodicidade semanal. Esse modelo foi computado por meio do pacote EpiEstim do R [10]. O código para reprodução dessa análise está disponível em rtlive_maringa.

Aqui, consideramos que o número de infectados diário tem a mesma interpretação do número de confirmados diário. Sabemos, porém, que existe um atraso na confirmação dos casos. Desde o dia em que a pessoa é infectada, leva-se um tempo para o vírus incubar (cerca de cinco dias) e os sintomas se manifestarem. Além disso, também pode existir um atraso de notificação e de espera pelos resultados dos testes. Esses problemas seriam resolvidos se tivéssemos os dados referentes ao atraso de notificação e de data do começo dos sintomas para cada caso confirmado, a fim de realizar medidas de nowcasting [12]. Além disso, o modelo não leva em consideração o volume de testagem variável (isto é, que o número de testes varia diariamente), o que pode causar uma subestimação do número de reprodução. Ainda assim, decidimos apresentar as estimativas do número de reprodução instantâneo. Diferentemente do número de reprodução efetivo, as estimativas são mais sensíveis às variações mais recentes dos números de casos confirmados diários. Sua interpretação deve ser aliada à análise de outros parâmetros aqui expostos.

Do ponto de vista da saúde pública, desejamos reduzir o número de reprodução efetivo/instantâneo $\mathcal{R}(t)$ (idealmente para um valor abaixo do limiar epidêmico) por meio de ações de contenção como distanciamento social e medidas de higiene. Como já dissemos, incorporar esses ingredientes em um modelo é uma tarefa praticamente infactível, já que não temos acesso ao grau de adoção dessas medidas pela população e nem mesmo sabemos estimar suas efetividades. Por isso, ao invés de fazermos previsões usando o modelo SIR, decidimos estimar a evolução de seus parâmetros e do número de reprodução efetivo/instantâneo à medida que os dados relacionados ao número de casos de COVID-19 vão sendo atualizados.

Modelo alternativo para o cálculo do número de reprodução efetivo

Recorremos ao método utilizado pelo website rt.live para o cálculo do número de reprodução efetivo $\mathcal{R}_{ef}(t)$. Podemos dizer que essa abordagem estima a curva do número de reprodução com maior probabilidade de gerar o número de infectados diários considerando toda a história da epidemia. Inicialmente, consideramos que o número de reprodução $\mathcal{R}_{ef}(t)$ caminha aleatoriamente, isto é, $$ \log \mathcal{R}_{ef}(t) \sim \mathcal{N}[R_{ef}(t-1), \sigma=0.1]\ , $$ em que $\sigma$ indica o quão rapidamente o número de reprodução pode variar diariamente. Além disso, supomos que existe uma parcela de indivíduos que iniciam a epidemia e estão distribuídos como $$ \text{Primeiros infectados} \sim Exp(\lambda=50)\ . $$ A partir disso, podemos simular a epidemia por meio da relação $$ I(t) = \mathcal{R}_{ef}(t)\sum_s I(t-s)w(s)\ , $$ sendo $I(t)$ o número de infectados no tempo $t$, $w(s)$ a distribuição de probabilidade do intervalo de tempo de geração e $s$ o número de dias anteriores. A distribuição de probabilidade $w(s)$ é a distribuição da diferença das datas de infecção entre casos primários e secundários. Utilizamos as estimativas de Nishiura et al. [7] para gerar a distribuição $w(s)$ com média $\mu=4.7$ dias e desvio padrão $\sigma=2.9$ dias, como utilizados no rt.live. Para incorporar o atraso da notificação e, além disso, o tempo de incubação do vírus (= 5 dias [8]) utilizamos informações do banco de dados SRAG 2020 disponibilizados pelo DATASUS. A Síndrome Respiratória Aguda Grave (SRAG) é um conjunto de sintomas (febre, tosse/dor de garganta, dispineia/saturação de oxigênio/desconforto respiratório) e tem notificação obrigatória no Sistema de Informação da Vigilância Epidemiológia da Gripe (SIVEP-Gripe). Esse sistema existe desde 2009 e foi criado com intuito de coletar informações sobre SRAG na epidemia de H1N1. O preenchimento é realizado através de uma ficha cadastral que inclui dados como data dos primeiros sintomas e data de notificação. Em 21/03/2020, o SARS-CoV-2 passou a integrar o painel de vírus testados para casos de SRAG. Por isso, utilizamos os dados de casos positivados de COVID-19 em Maringá nesse banco de dados para estimar o atraso entre os primeiros sintomas e a notificação do caso. O dado que possuímos é parcial, pois nem todos os casos de COVID-19 evoluem para SRAG, mas consideramos que é uma aproximação razoável. Uma melhor solução seria utilizar a base de dados integral para todos os casos confirmados na cidade se estes fossem disponibilizados pela Secretaria de Saúde, porém, este não é o caso. A partir do banco de dados SRAG 2020 e do tempo de incubação do vírus, estimamos uma distribuição de probabilidade que denominaremos $h(s)$ sendo $s$ o número de dias anteriores à confirmação do caso, incluindo também o tempo médio para obtenção do resultado do teste da COVID-19 em Maringá (cerca de dois dias). A partir dessa distribuição, realizamos uma convolução para estimar a data de infecção dos casos, ou seja, $$ I(t) = I_{c}(t)*h(s)\ , $$ em que $I_{c}(t)$ o número de casos confirmados no tempo $t$. Dessa forma, incorporamos o atraso de notificação e o tempo de incubação ao modelo. Para considerar o volume de testagem variável, podemos utilizar uma distribuição binomial negativa com um parâmetro de exposição para cada dia. O parâmetro de exposição leva em consideração as diferentes taxas de testagem diárias. Em outras palavras, ajustamos o número de casos pela quantidade de “oportunidade” que um evento possui, ou seja, $$ \text{contagem} = \text{taxa}\cdot\text{exposição}\ . $$ No caso da COVID-19, a exposição é a quantidade de testes diários. Se há mais testes, é claro que possivelmente o número de casos confirmados deva ser maior. Incorporando essa informação ao modelo temos que o número de infectados pode ser escrito como $$ I(t) \sim NegBin(\mu, \alpha)\ , $$ $$ \mu = I(t) \cdot \text{exposição}\ , $$ $$ \alpha \sim Gamma(\mu=6,\sigma=1)\ , $$ de acordo com esta definição da distribuição binomial negativa. É importante ressaltar que a inclusão de novos dados diários mudam as estimativas de datas anteriores, por isso, a cada dia teremos uma curva diferente (e mais precisa). Amostramos esse modelo utilizando o pacote PyMC3 do Python. O código para reprodução dessa análise está disponível em rtlive_maringa.

Referências

[1] S. Picoli Jr., J. J. V. Teixeira, H. V. Ribeiro, L. C. Malacarne, R. P. B. Santos, R.S. Mendes, Spreading Patterns of the Influenza A (H1N1) Pandemic. PLoS One 8, e17823 (2011).

[2] K. Wu, D. Darcet, Q. Wang, D. Sornette, Generalized logistic growth modeling of the COVID-19 outbreak in 29 provinces in China and in the rest of the world. arXiv: 2003.10376 (2020).

[3] G. Dattoli, E. Di Palma, S. Licciardi, E. Sabia, On the Evolution of Covid-19 in Italy: a Follow up Note. arXiv: 2003.12667 (2020).

[4] Y. Liu, A. A. Gayle, A. Wilder-Smith, J. Rocklv, The reproductive number of COVID-19 is higher compared to SARS coronavirus. Journal of Travel Medicine 27, taaa021 (2020). J. T. Wu, K. Leung, G. M. Leung, Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study. The Lancet 395, 689 (2020).

[5] B. J. Coburn, B. G. Wagner, S. Blower, Modeling influenza epidemics and pandemics: insights into the future of swine flu (H1N1). BMC Medicine 7 (2009).

[6] A. J. Stier, M. G. Berman, L. M. A. Bettencourt, COVID-19 attack rate increases with city size. arXiv: 2003.10376 (2020).

[7] H. Nishiura, N.M. Linton, A.R. Akhmetzhanov, Serial interval of novel coronavirus (COVID-19) infections. International Journal of Infectious Diseases 93, 284-286 (2020).

[8] S. A. Lauer, K. H. Grantz, Q. Bi, F. K. Jones, Q. Zheng, H. R. Meredith, A. S. Azman, N. G. Reich, J. Lessler, The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: Estimation and application. Annals of Internal Medicine 172, 577-582 (2020).

[9] C. Fraser, Estimating Individual and Household Reproduction Numbers in an Emerging Epidemic. PLOS ONE 2(8), e758 (2007).

[10] A. Cori, N. M. Ferguson, C. Fraser, S. Cauchemez, A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics. American Journal of Epidemiology 178, 1505–1512 (2013).

[11] S. Abbott, J. Hellewell, J. Munday, S. Funk, The transmissibility of novel Coronavirus in the early stages of the 2019-20 outbreak in Wuhan: Exploring initial point-source exposure sizes and durations using scenario analysis. Wellcome Open Research 5, 17 (2020).

[12] L. S. Bastos, T. Economou, M. F. C. Gomes, D. A. M. Villela, F. C. Coelho, O. G. Cruz, O. S. T. Bailey, C. T. Codeço, A modelling approach for correcting reporting delays in disease surveillance data. Statistics in Medicine 38, 4363-4377 (2019).