Model Description
*****************

Basic model
===========

We model time :math:`t` as an *integer* with the meaning of "days since day 0",
where day 0 can be chosen at will.

The primary dynamic variable :math:`\delta N_{\mathrm{inf}}(t)` is the number of people
being infected on day :math:`t`. The symbol :math:`\delta` indicates that this
is an increment per day, not the total number of cases. The total number of
people that were infected between times :math:`t_0` and :math:`t_1` is:

.. math::

   N_{\mathrm{inf}}(t_0,t_1) = \sum_{t=t_0}^{t_1} \delta N_{\mathrm{inf}}(t)

For a population without immunity, the mean number of people infected by a
single person is, by definition, given by the symbol :math:`R_0` (the basic
reproduction number). For a non-zero number of immune people :math:`N_{\mathrm{imm}}`,
the actual mean number of people infected by a single person :math:`R_e` equals

.. math::

   R_e(t) = R_0\; \left(1-\frac{N_{\mathrm{imm}}(t)}{N_{\mathrm{pop}}}\right)

where :math:`N_{\mathrm{pop}}` is the total number of people in the population.

The *mean* time between successive infections is called the *generation time*
:math:`t_g`. It means that, on average :math:`t_g` days after being infected, a
person infects :math:`R_e` further individuals with the virus. In practice this
"passing on the virus" is spread over a time period, with different
probabilities for different times :math:`t` after the day of infection
:math:`t_i`. We describe this "time spreading function of reproduction" as
:math:`P_r(t-t_i)`, with obeys the following normalization:

.. math::

   \sum_{t_i<t} P_r(t-t_i) = 1

because :math:`P_r(t-t_i)` is a probability density function (PDF).
The mean generation time is then its first moment:

.. math::

   t_g = \sum_{t_i<t} (t-t_i)\;P_r(t-t_i)

With :math:`P_r(t-t_i)` and :math:`R_e(t)` we can now construct the
*infectivity function* :math:`r(t_i,t)`

.. math::

   r(t_i,t) = R_e(t-1) \; P_r(t-t_i)

which is the mean number of persons being infected on day :math:`t` by an
individual who was infected on day :math:`t_i<t`. The :math:`t-1` in the
:math:`R_e(t-1)` is because :math:`R_e` is not yet computed for time
:math:`t` when it is used in the computation of :math:`r(t_i,t)`.

Note that the use of :math:`R_e(t-1)` for the computation of the infectivity
function is the only place where a non-linear feedback is involved. 

In the basic version of our model we make the following simplified assumption: A
person infected on day :math:`t_i` is assumed to be able to infect other people
during a period of days from :math:`g_0` up to (and including) :math:`g_1` days
after day :math:`t_i`. During that period, in this simple model, the infectivity
:math:`r` is constant, while before and after, the infectivity :math:`r` is
zero:

.. math::

   P_r(t-t_i) = \left\{\begin{matrix}
   \frac{1}{g_1-g_0+1} & \mathrm{if} \; g_0 \le t-t_i \le g_1 \\
    0 & \mathrm{otherwise}
   \end{matrix}\right.

and the generation time is simply :math:`t_g = (g_0+g_1)/2`. 

With :math:`r(t_i,t)` defined, the *epidemic dynamic equation* for :math:`\delta
N_{\mathrm{inf}}(t)` becomes

.. math::

   \delta N_{\mathrm{inf}}(t) = \sum_{t_i<t}\; r(t_i,t)\; \delta N_{\mathrm{inf}}(t_i)

The system of equations can be closed by assuming, for instance, that all (and
only those) people that have been infected are immune:

.. math::

   N_{\mathrm{imm}}(t) = N_{\mathrm{inf}}(-\infty,t)

and by assuming :math:`N_{\mathrm{pop}}` to be constant and given, e.g.
:math:`N_{\mathrm{pop}}=83,000,000` for Germany.

The integration of this epidemic dynamic equation is straightforward: time-stepping in
steps of one day. For :math:`R_0>1` this leads to an exponential growth of
:math:`\delta N_{\mathrm{inf}}(t)` that saturates as a result of herd immunity.

The dynamics of the epidemic is entirely described by :math:`\delta
N_{\mathrm{inf}}(t)`. Other quantities, such as the current number of ill
people, the number of people in the hospital and on an intensive care unit
(ICU), the number of people that have recovered or died, can all be computed
from :math:`\delta N_{\mathrm{inf}}(t)`.


Derived quantities
==================

For practical purposes, an epidemic dynamic model should be able to predict how many
people will, at some time :math:`t` in the future, be in a particular "state"
(e.g. having symptoms, being in the hospital, occupying an ICU bed, or having
deceased). Averaged over many different progressions of the disease, one can
define a probability function :math:`p_{\mathrm{state}}(t-t_i)` for a person
that has become infected on day :math:`t_i` to be in that state on day
:math:`t`. These probability functions (PF) are, however, *not normalized*,
because they are not probability density functions (PDFs), and because not every
individual will ever acquire such a state. However, by their nature as PFs, they
must obey

.. math::

   0\le p_{\mathrm{state}}(t-t_i)\le 1

at all :math:`t\ge t_i`, and :math:`p_{\mathrm{state}}(t-t_i)=0` for :math:`t<t_i`.
Typically the numbers one gets from databases or reports is the fraction of
people who acquire a certain state, e.g. that about 1% of reported cases get
admitted to an ICU. That information is not sufficient to estimate the number of
required ICU beds. The mean duration :math:`T_{\mathrm{ICU}}` of the stay at the
ICU, the mean time after infection :math:`t_{\mathrm{ICU}}` that person is
admitted to an ICU, and the probability :math:`p_{\mathrm{ICU}}` that the
person will enter such a state are required to make the simplest possible model
of :math:`p_{\mathrm{ICU}}(t-t_i)`. In general this simplest model would look as
follows:

.. math::

   p_{\mathrm{state}}(t-t_i) = \left\{\begin{matrix}
    p_{\mathrm{state}} & \mathrm{if} \; t_{\mathrm{state}} \le t-t_i \le t_{\mathrm{state}}+T_{\mathrm{state}} \\
    0 & \mathrm{otherwise}
   \end{matrix}\right.

where "state" can be replaced by "ICU", "ill", "hospital", "recovered" or any
other state of interest. More realistic expressions for
:math:`p_{\mathrm{state}}(t-t_i)` can be generated, for instance, by convolving
this function with a Gaussian to account for the variance in
:math:`t_{\mathrm{ICU}}` and :math:`T_{\mathrm{ICU}}`.

The computation of the number of people in that state at time :math:`t` is
then straightfoward:

.. math::

   N_{\mathrm{state}}(t) = \sum_{t_i<t}\; P_{\mathrm{state}}(t-t_i)\; \delta N_{\mathrm{inf}}(t_i)


Adding features to the basic model
==================================

It is relatively straightforward to add features to the basic model. Here we give a few
examples.


Effect of sudden policy changes
-------------------------------

The implementation of a (partial) lockdown on day :math:`t_{\mathrm{lock}}` will suddenly
change the :math:`R_0`. In EpiDemo this can easily be implemented. One asks the
model to integrate up to :math:`t_{\mathrm{lock}}`. Then one changes the :math:`R_0`
variable, and asks the model to continue integrating until the end.


Effect of a gradually implemented vaccination campaign
------------------------------------------------------

Already built into EpiDemo is the possibility to prescribe an array of values
for :math:`N_{\mathrm{vacc}}(t)` for the entire modeling time domain, which
represents the number of people that have been vaccinated up to that
point. EpiDemo then calculates :math:`N_{\mathrm{imm}}(t)` as follows:

.. math::

   N_{\mathrm{imm}}(t) = \min\left(N_{\mathrm{inf}}(-\infty,t)+N_{\mathrm{vacc}}(t),N_{\mathrm{pop}}\right)


Effect of waning immunity
-------------------------

There are several ways by which a waning natural immunity can be implemented.
The simplest being an exponential decay: Instead of writing
:math:`N_{\mathrm{imm}}(t) = N_{\mathrm{inf}}(-\infty,t)` one can write:

.. math::

   N_{\mathrm{imm}}(t) = \sum_{t'=-\infty}^{t} \delta N_{\mathrm{inf}}(t')\;e^{-(t-t')/\tau}

with :math:`\tau` being the immunity waning time scale. 


Multi-group version of the model
================================

A real society is heterogeneous. Certain communities stay mostly among each
other, having little interaction with other communities. An infection within
such a community will, in the ideal (and unrealistic) case of zero interaction
with the outside world, develop its own dynamics, with its own :math:`R_0`
value. Since the dynamics is exponential, even a small difference with the
value of :math:`R_0` outside that group will, after some time, lead to vast
differences in the fraction of infected people compared to the outside world.

Suppose a group is so isolated that it managed to avoid infection altogether,
even though already 5% of the population outside that group has become
infected. Suppose that the behavior in this group is such that contacts are much
more close and common than in the population at large (perhaps due to a false
sense of security). If, for whatever reason, a single individual within this
group *does* get infected, the behavior of this group may set their value of
:math:`R_0` to be much larger than that of the outside world. As a result, the
virus spreads like wildfire inside this group, and soon the fraction of infected
within the group greatly exceeds that of the more careful population outside.

If that group is large enough, the total number of registered cases in the
country may suddenly see a sharp rise seemingly coming out of nowhere, as the
total reported new infections of the group briefly outperforms the rest of the
country, before the group reaches herd immunity or changes its behavior.

Since a group is never perfectly isolated, once the number of cases in the group
massively dominates that of the outside world, even a small level of cross
infection to the outside can then dominate the outside infection rates as
well. 

In Germany the `massive outbreak occurred at a slaughterhouse in the region of
Gütersloh <https://de.m.wikipedia.org/wiki/Tönnies_Holding#Covid_Stammwerk>`_ in
the summer of 2020 is likely an example of this, where the virus could, over an
extended period of time, multiply in a relatively isolated community until a
large fraction of that community was infected (see e.g. `here
<https://www.faz.net/aktuell/gesellschaft/gesundheit/coronavirus/expertin-ausbruch-bei-toennies-schwelt-sehr-wahrscheinlich-schon-laenger-16820726.html>`_).
In the RKI data it can be clearly seen as a "bump" in the total number of
reported new infections around June, which then spread to other regions of
Germany and dissipated as a result of the overall negative :math:`R_0-1` in
Germany at that time.

This shows that even if one group is relatively isolated from the rest of
society, if an epidemic is out of control in this group, the epidemic dynamics of
society can become "enslaved" by the epidemic dynamics within the semi-isolated
group (see e.g. the Keeling and Rohani book).

Note that this is different from a "superspreading event", in which a single (or
a few) people infect a very large number of people all at once due to
circumstances that are exceptionally conducive to viral spread. 

Most "groups" within a society are not so strongly isolated, and if contacts
between groups are very common, there is no need to treat the groups separately
in a model. The mean :math:`R_0` value will then be enough to make reliable
predictions.

But there are intermediate cases that are of particular interest. For instance,
the attempts to isolate the elderly. Particularly the elderly in nursing homes
are, out of necessity, in regular contact with nursing personnel. The degree of
care taken to avoid cross infection from nursing personnel to patients
determines the rate of cross-infection. How careful do they have to be to be
successful in case the epidemic gets out of hand? Another example is the role
of subcommunities of young people, some of whom may behave less careful than
the rest of the population. How dangerous is this? None of these questions
can be easily quantitatively answered, because the knowledge of the parameters
is unavailable. But some simple scenarios can be tested.

In the multi-group model, :math:`\delta N_{\mathrm{inf}}(t)` is replaced by
:math:`\delta N_{\mathrm{inf},k}(t)`, where :math:`k` is the index of the
group. Python counts array indices from 0, so :math:`k\;\in [0,n-1]` with
:math:`n` being the number of groups. The epidemic dynamic equation then becomes

.. math::

   \delta N_{\mathrm{inf},k}(t) = \sum_{l=0}^{n-1} \sum_{t_i<t}\;
   r_{kl}(t_i,t)\; \delta N_{\mathrm{inf},l}(t_i)

where :math:`r_{kl}(t_i,t)` is the *infectivity function* for infection from
group :math:`l` to group :math:`k`:

.. math::

   r_{kl}(t_i,t) = R_{\mathrm{e},kl}(t) \; P_r(t-t_i)

where we assume that :math:`P_r(t-t_i)` is the same for all groups. Note the
order of the indices: :math:`k\leftarrow l`, which is chosen to conform to the
standard matrix multiplication index order. The "reproduction matrix" can
be computed from the "basic reproduction matrix" :math:`R_{0,kl}`

.. math::

   R_{\mathrm{e},kl}(t) = R_{0,kl}\; \left(1-\frac{N_{\mathrm{imm},k}(t)}{N_{\mathrm{pop},k}}\right)

The values of the matrix :math:`R_{0,kl}` have to be specified. The sums
of the columns

.. math::

   R_{0,l} = \sum_{k=0}^{n-1} R_{0,kl}

give the the basic reproduction number for the average individual within group
:math:`l`, where :math:`R_{0,ll}` is the infection within the group and
:math:`R_{0,l}-R_{0,ll}` the infection to the other groups. 

The off-diagonal elements of :math:`R_{0,kl}` are tricky, because they depend
not only on the behavior of people within the group, but also on how many people
from the other groups they interact with. That depends also on the sizes of the
groups. The matrix is thus also not necessarily symmetric. If an individual from
group :math:`k` meets with an individual from group :math:`l`, and both equally
much expose the other, then :math:`R_{0,kl}N_{\mathrm{pop},l}
=R_{0,lk}N_{\mathrm{pop},k}`. 

For the rest, the model functions exactly the same as the one-group model.