# Epidemics: towards understanding
undulation and decay

Niko Sauer
University of Pretoria, South Africa

###### Abstract

Undulation of infection levels, usually called waves, are not well understood.
In this paper we propose a mathematical model that exhibits undulation and
decay towards a stable state. The model is a re-interpretation of the original
SIR-model obtained by postulating different constitutive relations whereby
classical logistic growth with recovery is obtained. The recovery relation is
based on the premise that infectiousness only lasts for some time. This leads
to a differential-difference (delay) equation which intrinsically exhibits
periodicity in its solutions but not necessarily decay to asymptotically
stable equilibrium. Limit cycles can indeed occur. An appropriate
linearization of the governing equation provides a firm basis for heuristic
reasoning as well as confidence in numerical calculations.

> MSC: 92D30 — Epidemiology; 34K13 — Periodic solutions.

> For my soul-mate Adri Prinsloo (1974–2021) whose penetrating questions —why,
> not how— contributed immensely to the development and understanding of this
> work. She left this life too soon.

##  1 Introduction

In a pioneering paper Hutchinson [10] states: “…that circular paths often
exist which tend to be self-correcting within certain limits, but which break
down, producing violent oscillations …” An equation to model this situation,
is given in a footnote to the paper as

$\frac{dy}{dt}(t)=y(t)[1-y(t-\tau)]$ (1.1)

with $\tau$ a ‘time lag’. This, with some constants added, came to be known as
Hutchinson’s equation. The oscillatory nature of its solutions was the subject
of a number of mathematical studies, the earliest of which are Cunningham [5],
Wright [18] and Jones [11], sometimes in an equivalent form. A generalization
dealing with several time lags is treated in Gopalsamy [9]. This was anteceded
by Van der Plank [15], [16] who considered plant diseases in which dormant as
well as an infectious periods are taken into account. The equation that
carries his name, suitably transformed, is

$\frac{dy}{dt}(t)=Ry(t)[y(t-\tau_{{}_{{}^{d}}})-y(t-\tau_{{}_{{}^{i}}})],$
(1.2)

which is of interest to us when the period $\tau_{{}_{{}^{d}}}$ of dormancy is
taken to be zero. By and large, studies of equations such as these have
concentrated on long term behaviour of solutions, particularly decay to a
point of equilibrium; “flattening of the curve”. A phenomenological model for
the prediction of “waves” is presented by Cacciapaglia, Cot & Sannino in [4].

In present-day clinical contexts, the phenomenon of “waves of infection” seems
to be very much at the forefront but does not appear to be well-understood or
even defined. It is the purpose of this paper to align this with Hutchinson’s
“violent oscillations” which have been mathematically shown to be exhibited by
equations such as (1.1) and its generalizations. The term wave may be
inappropriate since wave phenomena inseparably involve both time and space.
For this reason we have chosen the word undulation. This phenomenon is known
to occur in plant as well as animal populations.

Instead of (1.1) we shall use the “logistic delay equation” as stated in Ruan
[14]. It is of the form
$x^{\prime}(t)=rx(t)[1-a_{{}_{{}^{1}}}x(t)-a_{{}_{{}^{2}}}x(t-\tau)]$, with
prime denoting the time derivative $d/dt$. This equation turns up occasionally
without indication of the assumptions made to derive it. In §2 we give a
systematic derivation based on a general view of the SIR model introduced by
Kermack & McKendick [12]. This view also leads to the “theta-model” which can
be used to obtain better correspondence to observed data. Moreover, this
approach establishes parameters with specific significance which is easily
lost when the treatment is entirely mathematical. Section 3 deals with
normalization and scaling of the equations to obtain more familiar forms. It
is also shown there that the theta-model equation can be transformed to the
‘standard’ form with constants and variables having different meanings.

Some general results are obtained in §4 such as positivity of solutions when
the ‘initial history’ is so, and an upper bound which implies that solutions
cannot grow in an unbridled manner. In §5 we obtain a ‘natural’ linearization
and in §6 treat the linear homogeneous problem with the aid of the Laplace
transform. This elaborate treatment serves to augment rather sketchy
treatments found in the literature, constantly keeping track of the model
parameters. Inversion of the Laplace transform is discussed in §7. This leads
to a series representation of the solution of the linear homogeneous problem.
Section 8 gives sharp estimates of the position of poles and terms in the
series solution. It is shown that convergence to equilibrium is guaranteed if
a constant, given in terms of the parameters, is sufficiently small.

The significance of the preceding analysis for the nonlinear problem is
discussed in §9. In §10 we present a numerical example to illustrate that
undulations, as exhibited by the model, can qualitatively be in accordance
with observed phenomena. The example shows that decay to an equilibrium point
is possible, but also that limit cycles are possible steady states. The
numerical procedure, based on a construction in §4, is given in §11.

The concluding remarks of §12, although unscientific, have some seriousness
about them.

##  2 SIR, Verhulst and more

Fundamental to many mathematical descriptions of epidemics is the SIR model
and its variants. Here three quantities, the number of susceptibles $S(t)$,
the number of infectious individuals $I(t)$, and the number recovered
(restored?) $R(t)$ at time $t$, are related by the ‘conservation principle’

$S(t)+I(t)+R(t)=N$ (2.1)

with the constant $N$ denoting the total ‘population’ considered. The dynamics
of the epidemic is in the system of ordinary differential equations

$\displaystyle I^{\prime}(t)=-[S^{\prime}(t)+R^{\prime}(t)];$ (2.2)
$\displaystyle S^{\prime}(t)=\mathfrak{F}_{{}_{{}^{S}}}(S(t),I(t),R(t));$
(2.3) $\displaystyle
R^{\prime}(t)=\mathfrak{F}_{{}_{{}^{R}}}(S(t),I(t),R(t)).$ (2.4)

Equation (2.2) is simply (2.1) differentiated. For specific purposes the
‘driving forces’ $\mathfrak{F}_{{}_{{}^{S}}}$ and $\mathfrak{F}_{{}_{{}^{R}}}$
are chosen by postulates known as constitutive relations. In the original SIR-
model, the postulates are: $\mathfrak{F}_{{}_{{}^{S}}}=-\beta I(t)[S(t)/N]$
and $\mathfrak{F}_{{}_{{}^{R}}}=\gamma I(t)$ with $\beta$ and $\gamma$
positive ‘rate’ constants (dimension [time]${}^{{}^{{}_{-1}}}$).

The Verhulst logistic growth model [17], originally aimed at population
growth, is sometimes used for epidemics when recoveries do not occur. It may
be considered as a fundamental principle of population dynamics. In this model
the population size $N$ is called the carrying capacity, the number of
individuals that can be infected, and recovery is ignored. Let $P(t)$ be the
probability of finding an infectious individual at time $t$. Then the relevant
constitutive relation, to which is added an assignment of probability, is
expressed as follows:

$\displaystyle\mathfrak{F}_{{}_{{}^{S}}}=-\beta I(t)[1-P(t)];$ (2.5)
$\displaystyle P(t)=I(t)/N,$ (2.6)

where $1-P(t)$ signifies the probability of finding a ‘healthy’ individual.
This, combined with (2.2), yields the classical logistic equation

$I^{\prime}(t)=\beta I(t)[1-I(t)/N]$

which is at the core of many informed speculations.

We now turn to the situation where recovery is also taken into account and
postulate the relation

$\mathfrak{F}_{{}_{{}^{R}}}=\gamma I(t)P(t-\tau),$ (2.7)

with $\tau>0$ the period of infectiousness. This means that the rate of
recovery is proportional to the probability of infectiousness occurring at the
earlier instant $t-\tau$. If $P(t)$ is still specified according to (2.6), the
relations (2.5), (2.7) gives the logistic-recovery equation

$I^{\prime}(t)=I(t)\\{\beta[1-(I(t)/N)]-\gamma(I(t-\tau)/N)\\}.$ (2.8)

This is the equation we shall study, although expressed differently.

In a model for competing species Gilpin and Ayala [8] essentially chose
$P(t)=[I(t)/N]^{\theta}$ with $\theta$ a positive constant. In the context of
our discussion here, this gives rise to the theta-logistic-recovery equation

$I^{\prime}(t)=I(t)\\{\beta[1-(I(t)/N)^{\theta}]-\gamma(I(t-\tau)/N)^{\theta}\\}.$
(2.9)

One may view this as a mathematical generalization of the logistic model
($\theta=1$) to manipulate the sigmoidal curve, but it can be grounded in
probability theory (Feller [7, II.5–8, Randomization]).

Constitutive relations in SIR models can be found in Della Morte & Sannino [6]
and Buonomo & Cerasuolo [3]. In the latter a time lag is introduced in
$\mathfrak{F}_{{}_{{}^{S}}}$. A SIR-model with delay in
$\mathfrak{F}_{{}_{{}^{R}}}$, not the same as (2.7), has been suggested by
Reiser [13].

##  3 Normalization, scaling and a reduction

We first normalize the equation (2.8) by the setting $F(t)=I(t)/N$ so that
$I=N$ corresponds to $F=1$, and $P(t)=F(t)$. One may think of $F(t)$ as the
level of infectiousness at time $t$. The result is

$F^{\prime}(t)=F(t)\\{\beta[1-F(t)]-\gamma F(t-\tau)\\}.$ (3.1)

This is a differential-difference equation or delay equation. To solve the
equation for times $t>0$ we need to know the state (history) of $F$ for times
$t\in[-\tau,0]$ (see e.g., Bellman & Cooke [2]). Thus we have the initial
condition

$F(t)=\Phi(t)\ \text{for }-\tau\leq t\leq 0,$ (3.2)

with $\Phi$ a given function defined on $[-\tau,0]$.

If it is assumed that an asymptotically, nonzero stable state for solutions of
(3.1) exists, that is, if $F^{\prime}(t)\to 0$ as $t\to\infty$ and the limit
$F_{{}_{{}^{\infty}}}=\lim_{t\to\infty}F(t)$ exists, it follows that
$F_{{}_{{}^{\infty}}}=\beta/(\beta+\gamma)$. We use this parameter to scale
the equation (3.1). In addition the parameter $\tau$ is used as unit of time
to obtain a completely dimensionless formulation. Thus we define the new
variables $t_{{}_{{}^{\tau}}}=t/\tau$ and
$f(t_{{}_{{}^{\tau}}})=F(t)/F_{{}_{{}^{\infty}}}=F(\tau
t_{{}_{{}^{\tau}}})/F_{{}_{{}^{\infty}}}$ to obtain, in place of (3.1),

$\displaystyle f^{\prime}(t_{{}_{{}^{\tau}}})$ $\displaystyle=\beta\tau
f(t_{{}_{{}^{\tau}}})\left\\{\left[1-\left(\frac{\beta}{\beta+\gamma}\right)f(t_{{}_{{}^{\tau}}})\right]-\left(\frac{\gamma}{\beta+\gamma}\right)f(t_{{}_{{}^{\tau}}}-1)\right\\}$
$\displaystyle=\tilde{\beta}f(t_{{}_{{}^{\tau}}})\\{[1-\beta^{{}^{{}_{*}}}f(t_{{}_{{}^{\tau}}})]-\gamma^{{}_{{}^{*}}}f(t_{{}_{{}^{\tau}}}-1)\\}$
(3.3)

with $\tilde{\beta}:=\beta\tau$, $\beta^{{}^{{}_{*}}}:=\beta/(\beta+\gamma)$,
$\gamma^{{}_{{}^{*}}}:=\gamma/(\beta+\gamma)$ and the prime denoting
$d/dt_{{}_{{}^{\tau}}}$. We note also the identity
$\beta^{{}^{{}_{*}}}+\gamma^{{}_{{}^{*}}}=1$ which will be important at a
later stage and gives rise to an alternative form discussed in §12. The
initial condition (3.2) now has the form

$f(t_{{}_{{}^{\tau}}})=\phi(t_{{}_{{}^{\tau}}}):=\frac{\Phi(\tau
t_{{}_{{}^{\tau}}})}{F_{{}_{{}^{\infty}}}}\ \text{for }-1\leq
t_{{}_{{}^{\tau}}}\leq 0.$ (3.4)

It is clear that if a non-zero stable asymptote exists, it has the value
$f_{{}_{{}^{\infty}}}=1$. Also note that (3.3) has the same form as the
“logistic delay equation” mentioned in §1 but the constants can be interpreted
in terms of the constitutive relations (2.5), (2.7).

Under the transformations just described, the equation (2.9) has the form

$f^{\prime}(t_{{}_{{}^{\tau}}})=\tilde{\beta}f(t_{{}_{{}^{\tau}}})\\{1-\beta^{{}^{{}_{*}}}f^{\theta}(t_{{}_{{}^{\tau}}})-\gamma^{{}_{{}^{*}}}f^{\theta}(t_{{}_{{}^{\tau}}}-1)\\}.$

This equation can assume a less intimidating look if we replace
$f^{\theta}(t_{{}_{{}^{\tau}}})$ with $f(t_{{}_{{}^{\tau}}})$ to obtain an
equation of precisely the same form as (3.3) except that now
$\tilde{\beta}=\theta\beta\tau$. Thus we go on to study (3.3).

##  4 Some general considerations

In this section we consider the scaled equation (3.3) together with the
initial condition (3.4). We shall abuse notation by writing $t$ instead of the
dimensionless time $t_{{}_{{}^{\tau}}}$. The initial level
$f_{{}_{{}^{0}}}:=f(0)$ will be of significance.

First we obtain a formal representation of the solution by the substitution
$g(t)=1/f(t)$, familiar for equations of the Bernoulli-kind. This leads to

$g^{\prime}(t)+\tilde{\beta}[1-\gamma^{{}_{{}^{*}}}f(t-1)]g(t)=\beta^{{}^{{}_{*}}}\tilde{\beta}.$
(4.1)

If we interpret the term $f(t-1)$ as $1/g(t-1)$ this is a differential-
difference equation for $g$ with initial condition

$g(t)=1/\phi(t)\text{ for }t\in[-1,0].$ (4.2)

Associated with this equation we have the integrating factor

$i(t):=\tilde{\beta}[t-\gamma^{{}_{{}^{*}}}\psi(t)];\ \text{ with
}\psi(t):=\int_{0}^{t}f(s-1)\kern 1.00006ptds=\int_{-1}^{t-1}f(s)\kern
1.00006ptds,$ (4.3)

so that

$\tfrac{d}{dt}[\exp\\{i(t)\\}g(t)\\}]=\beta^{{}^{{}_{*}}}\tilde{\beta}\exp\\{i(t)\\}.$
(4.4)

This can be integrated directly. However, further integration by parts of the
term on the right yields:

$\displaystyle\exp\\{i(t)\\}g(t)$
$\displaystyle=g_{{}_{{}^{0}}}-\beta^{{}^{{}_{*}}}+\beta^{{}^{{}_{*}}}\exp\\{i(t)\\}$
$\displaystyle\qquad+\tilde{\beta}\beta^{{}^{{}_{*}}}\gamma^{{}^{{}_{*}}}\int_{0}^{t}\exp\\{i(s)\\}f(s-1)\kern
1.00006ptds,$ (4.5)

where $g_{{}_{{}^{0}}}=1/f_{{}_{{}^{0}}}$.

The formal calculations above can be placed on a firmer footing under the
following hypotheses about the initial state which will be taken for granted
from now on:

  1. H1.

The initial state $\phi$ is continuous on $[-1,0]$.

  2. H2.

$\phi(t)>0$ for $t\in[-1,0]$.

We immediately note that the integrating factor $i(t)$ exists and is
differentiable for $t>0$.

The initial value problem we study is well-posed in the sense of the following
result.

######  Theorem 4.1.

Under the assumptions H1 and H2:

  1. (a)

The function $g(t)$ as represented in (4.5) is a positive solution of (4.1),
(4.2) with $g(0)=g_{{}_{{}^{0}}}$.

  2. (b)

The function $f(t)=1/g(t)$ is a positive solution of (3.3), (3.4).

###### Proof.

The proof of the two assertions will simultaneously unfold by progression over
the time intervals $[0,1),[1,2),\dots$ as in many instances to be found in
[2].

We begin with $0\leq t<1$. Here the function $\psi$ is totally determined by
the initial state $\phi$; it is in fact differentiable and positive. Thus the
integrating factor is of suitable nature and the function $g(t)$, as
determined by (4.5), is positive and solves the initial value problem (4.1),
(4.2). It follows that $f(t)>0$ solves (3.3), (3.4). Thus (a) and (b) are
established in $[0,1)$. In addition the limit as $t\to 1$ defines the
functions $g(1)$ and $f(1)$

The same argument can be followed in the interval $[1,2)$, as the crucial
properties have been established in $[0,1]$. It is clear that a formal
induction argument will lead to the required outcome. ∎

The next result shows that there are limits to infectivity levels and that
there is at most one stable asymptote.

######  Theorem 4.2.

Under the assumptions H1 and H2 the solution of the initial value problem
(3.3), (3.4) is restricted in the following ways:

  1. (a)

It is bounded. Specifically,

$0<f(t)<\frac{f_{{}_{{}^{0}}}}{\beta^{{}^{{}_{*}}}f_{{}_{{}^{0}}}+(1-\beta^{{}^{{}_{*}}}f_{{}_{{}^{0}}})\exp\\{-\tilde{\beta}t\\}}\text{
for }t>0.$

  2. (b)

If the limit $f_{{}_{{}^{\infty}}}=\lim_{t\to\infty}f(t)$ exists, it is equal
to 1

###### Proof.

Since $f(t)>0$ for $t\geq-1$, we see from (4.3) that $i(t)<\tilde{\beta}t$.
From the representation (4.5) follows that

$g(t)>\beta^{{}^{{}_{*}}}+(g_{{}_{{}^{0}}}-\beta^{{}^{{}_{*}}})\exp\\{-i(t)\\}>\beta^{{}^{{}_{*}}}+(g_{{}_{{}^{0}}}-\beta^{{}^{{}_{*}}})\exp\\{-\tilde{\beta}t\\}.$

The stated upper bound is obtained by reciprocation and further manipulation.

To prove (b) let us assume that the limit is zero. Then for given
$\varepsilon\in(0,1)$ there exists $t_{{}_{{}^{0}}}$ such that for
$t>t_{{}_{{}^{0}}}$ both $f(t)$ and $f(t-1)$ are less than $\varepsilon$. This
implies that $f^{\prime}(t)>\tilde{\beta}f(t)[1-\varepsilon]$ since
$\beta^{{}^{{}_{*}}}+\gamma^{{}_{{}^{*}}}=1$. It now follows that
$f(t)>f(t_{{}_{{}^{0}}})\exp\\{\tilde{\beta}[1-\varepsilon](t-t_{{}_{{}^{0}}})\\}$
for all $t>t_{{}_{{}^{0}}}$ which is absurd.

It is seen from (3.3) that the limit $\lim_{t\to\infty}f^{\prime}(t)$ exists.
The mean value theorem shows that this limit is zero, and hence
$f_{{}_{{}^{\infty}}}[1-f_{{}_{{}^{\infty}}}]=0$. ∎

We note that $f(t)<M:=\max\\{f_{{}_{{}^{0}}},1/\beta^{{}^{{}_{*}}}\\}$.

##  5 Linearization

For the problem at hand we arrive at a suitable linearization by shifting the
(expected) equilibrium level from $f=1$ to $f=0$. To avoid an undue
proliferation of symbols, we once again abuse notation by the replacing
$f(t)-1$ by $f(t)$. The governing equation (3.3) then has the form

$f^{\prime}(t)+\tilde{\beta}[\beta^{{}^{{}_{*}}}f(t)+\gamma^{{}_{{}^{*}}}f(t-1)]=-\tilde{\beta}f(t)[\beta^{{}^{{}_{*}}}f(t)+\gamma^{{}_{{}^{*}}}f(t-1)].$
(5.1)

Of course, the initial condition (3.4) is adapted accordingly.

One immediately notes that the left of (5.1), as opposed to the right, is
linear. Intuitively, if $f(t)$ is close to zero for large $t$, the nonlinear
term which is quadratic in $f$ will be insignificant in the long run.
Formally, the right of (5.1) linearizes to zero. This leads us to the linear
homogeneous equation

$f^{\prime}(t)+\tilde{\beta}[\beta^{{}^{{}_{*}}}f(t)+\gamma^{{}_{{}^{*}}}f(t-1)]=0,$
(5.2)

which will be studied in detail.

Some simplifying notation is introduced:
$b:=\tilde{\beta}\beta^{{}^{{}_{*}}}$, $c:=\tilde{\beta}\gamma^{{}_{{}^{*}}}$
and $B:=c\exp\\{b\\}$.

We obtain an estimate for solutions of the non-homogeneous problem

$\left.\begin{aligned} &f^{\prime}(t)+bf(t)+cf(t-1)=v(t)\text{ for }t>0;\\\
&f(t)=\phi(t)\text{ for }-1\leq t\leq 0,\end{aligned}\right\\}$ (5.3)

with $v$ a given function, continuous on $[0,\infty)$.

######  Theorem 5.1.

Every solution of (5.3) satisfies

$|f(t)|\leq[F_{{}_{{}^{0}}}+V(t)]\exp\\{(B-b)t\\};\quad t\geq 0,$

with
$F_{{}_{{}^{0}}}:=|f_{{}_{{}^{0}}}|+c\int_{-1}^{0}\exp\\{bs\\}|\phi(s)|\kern
1.00006ptds$ and $V(t)=\int_{0}^{t}\exp\\{bs\\}|v(s)|\kern 1.00006ptds$.

###### Proof.

Let us write the first equation in (5.3) in the form

$\tfrac{d}{dt}[\exp\\{bt\\}f(t)]=\exp\\{bt\\}[v(t)-cf(t-1)],$

integration of which yields

$\displaystyle[\exp\\{bt\\}f(t)]=f_{{}_{{}^{0}}}-c\int_{0}^{t}\exp\\{bs\\}f(s-1)\kern
1.00006ptds+\int_{0}^{t}\exp\\{bs\\}v(s)\kern 1.00006ptds$
$\displaystyle\qquad=f_{{}_{{}^{0}}}-c\int_{-1}^{0}\exp\\{bs\\}\phi(s)\kern
1.00006ptds-c\exp\\{b\\}\int_{0}^{t-1}\exp\\{bs\\}f(s)\kern 1.00006ptds$
$\displaystyle\qquad\qquad+\int_{0}^{t}\exp\\{bs\\}v(s)\kern 1.00006ptds.$

With $F(t):=\exp\\{bt\\}|f(t)|$, we obtain the integral inequality

$F(t)\leq F_{{}_{{}^{0}}}+V(t)+B\int_{0}^{t}F(s)\kern 1.00006ptds.$

Since $V$ is monotonically increasing, the Grönwall-Bellman inequality [1],
[2, Lemma 3.1] applies, i.e., $F(t)\leq[F_{{}_{{}^{0}}}+V(t)]\exp\\{Bt\\}$,
and the proof is complete. ∎

This result shows that the solution of (5.3) is unique and depends
continuously on the initial data.

##  6 The linear homogeneous equation

Our attention now turns to the linear equation (5.2) under the initial
condition $f(t)=\phi(t)\text{ for }-1\leq t\leq 0;f(0)=f_{{}_{{}^{0}}}$. The
approach is by the Laplace transform defined as
$\hat{f}(s):=\int_{0}^{\infty}\exp\\{-st\\}f(t)\kern 1.00006ptdt$ for complex
$s$ with positive real part. In fact, Thm. 5.1 with $v\equiv 0$ shows that for
$\text{Re}\,s>B-b$ the Laplace transform of (5.2) may be taken. The result is

$\left.\begin{aligned} &h(s)\hat{f}(s)=H(s);\\\ &h(s):=s+b+c\exp\\{-s\\};\\\
&H(s):=f_{{}_{{}^{0}}}-c\exp\\{-s\\}\int_{-1}^{0}\exp\\{-st\\}\phi(t)\kern
1.00006ptdt.\end{aligned}\right\\}$ (6.1)

For the inversion of $\hat{f}$ we need to study the zeros of the complex-
valued function $h(s)$ for $s=x+iy$ ($i=\sqrt{-1}$) in the complex plane. For
this the real and imaginary parts of $h$ must vanish and we have the equations

$\displaystyle(x+b)+c\exp\\{-x\\}\cos y=0;$ (6.2) $\displaystyle
y-c\exp\\{-x\\}\sin y=0,$ (6.3)

both of which need to be satisfied.

As a first step we eliminate the trigonometric terms in (6.2), (6.3) to obtain

$y^{2}=c^{2}\exp\\{-2x\\}-(x+b)^{2}.$ (6.4)

This defines a curve on which the solution points $s=x+iy$ must lie, but not
every point on the curve is necessarily a solution of (6.2) and (6.3).
Moreover, the curve so obtained is only defined for those $x$ for which the
right of (6.4) is non-negative. We investigate this question first.

######  Lemma 6.1.

With $B$ a positive constant:

  1. (a)

The equation $B\exp\\{-u\\}=u$ has a unique positive solution
$u_{{}_{{}^{M}}}$ and $B/(1+B)\leq u_{{}_{{}^{M}}}<B$.

  2. (b)

If $B>e^{-1}$ the equation $B\exp\\{u\\}=u$ has no positive solution.

###### Proof.

Assertion (a): The existence, uniqueness, positivity and the upper bound of a
solution is straightforward. To obtain the lower bound we note that
$\exp\\{-u\\}\geq 1-u$ for $u>0$. Hence
$u_{{}_{{}^{M}}}=B\exp\\{-u_{{}_{{}^{M}}}\\}\geq B(1-u_{{}_{{}^{M}}})$ and the
result follows.

Assertion (b) follows from the inequality $u\exp\\{-u\\}\leq e^{-1}$. ∎

In accordance with (6.4), let
$k^{{}^{{}_{(2)}}}\kern-3.0pt(x):=c^{2}\exp\\{-2x\\}-(x+b)^{2}=\break[c\exp\\{b\\}]^{2}\exp\\{-2(x+b)\\}-(x+b)^{2}.$

######  Proposition 6.1.

Let $B=c\exp\\{b\\}$ and $u_{{}_{{}^{M}}}$ as in Lemma 6.1. Then

  1. (a)

For $-b\leq x\leq x_{{}_{{}^{M}}}:=u_{{}_{{}^{M}}}-b$,
$k^{{}^{{}_{(2)}}}\kern-3.0pt(x)\geq 0$ with equality only if
$x=x_{{}_{{}^{M}}}$.

  2. (b)

For $x+b<0$ the function $k^{{}^{{}_{(2)}}}\kern-3.0pt(x)>0$ if $B>e^{-1}$.

###### Proof.

For statement (a), let $u=x+b\geq 0$. Then
$k^{{}^{{}_{(2)}}}\kern-3.0pt(x)=B^{2}\exp\\{-2u\\}-u^{2}=[B\exp\\{-u\\}+u][B\exp\\{-u\\}-u]$
which, by Lemma 6.1(a) is non-negative for $u\leq u_{{}_{{}^{M}}}$. To prove
(b) let $u=-(x+b)>0$. Then
$k^{{}^{{}_{(2)}}}\kern-3.0pt(x)=B^{2}\exp\\{2u\\}-u^{2}$. The result follows
from Lemma 6.1(b). ∎

It is of importance to introduce the parameter
$\rho:=\gamma/\beta=\gamma^{{}_{{}^{*}}}/\beta^{{}^{{}_{*}}}$. It corresponds
to the reciprocal of the basic reproduction number
$R_{{}_{{}^{0}}}=\beta/\gamma$ in the classical SIR model. This leads to the
relation $c=\rho b$. From now on we make an assumption stronger than suggested
by Lemma 6.1(b) namely

$B=c\exp\\{b\\}=\rho b\exp\\{b\\}>1.$ (6.5)

Under this assumption the function $k^{{}^{{}_{(2)}}}\kern-3.0pt(x)\geq 0$ for
$x\leq x_{{}_{{}^{M}}}$ so that $k^{{}^{{}_{(2)}}}\kern-3.0pt(x)$ can actually
be considered a square. We define the function $k$ by

$\displaystyle k^{2}(x):=k^{{}^{{}_{(2)}}}\kern-3.0pt(x)$
$\displaystyle=c^{2}\exp\\{-2x\\}-(x+b)^{2}$
$\displaystyle\qquad=B^{2}\exp\\{-2(x+b)\\}-(x+b)^{2}$
$\displaystyle\qquad=\rho^{2}b^{2}\exp\\{-2x\\}-(x+b)^{2};\quad x\leq
x_{{}_{{}^{M}}}.$ (6.6)

The curve $K$ defined in the complex plane by $s=x\pm ik(x)$ will be our next
concern. Because of symmetry we shall deal mainly with the positive branch
which will also be referred to as $K$.

######  Proposition 6.2.

The curve $K$ has the following properties:

  1. (a)

$k(x)\to\infty$ as $x\to-\infty$.

  2. (b)

At $x<x_{{}_{{}^{M}}}$ the tangent is negative.

  3. (c)

$x_{{}_{{}^{M}}}$ is positive if and only if $\gamma>\beta$, i.e., $\rho>1$.

###### Proof.

The assertion (a) is, by (6.6), straightforward.

From (6.6) we see that $k(x)k^{\prime}(x)=-[B^{2}\exp\\{-2(x+b)\\}+(x+b)]$. If
$x+b>0$, the term in brackets is positive. Since $k(x)>0$ it follows that
$k^{\prime}(x)<0$ for such $x$. If $x+b<0$ let $u=-(x+b)$ and it follows that
$B^{2}\exp\\{-2(x+b)\\}+(x+b)=B^{2}\exp\\{2u\\}-u>B^{2}(1+2u)-u>1+u>0$, by the
assumption (6.5). Thus (b) is established.

To prove (c) we notice that $k^{{}^{{}_{(2)}}}\kern-3.0pt(0)=c^{2}-b^{2}$.
Hence if $c>b$, the point $(0,[c^{2}-b^{2}]^{1/2})$ is on the curve $K$. From
(b) we see that $x_{{}_{{}^{M}}}$, where $k(x)=0$, must be positive. But $c>b$
means the same as $\gamma>\beta$. This argument can also be reversed. ∎

From what we have established so far, the following is significant:

  1. A.

Under the assumption (6.5) the zeros of the function $h$ are all to the left
of the vertical line $x=B$.

  2. B.

If $\beta\geq\gamma$ the zeros in question are to the left of the line
$x=\delta$ for arbitrary $\delta>0$. In fact, if $\beta=\gamma$,
$x_{{}_{{}^{M}}}=0$.

Thus inversion of the Laplace transform becomes a distinct possibility if the
zeros on the curve $K$ can be located. For that we need to obtain information
about the points $(x,k(x))$ for which the equations (6.2) and (6.3) are
actually satisfied. Towards this we eliminate the exponential terms from these
equations to obtain the relation $y+(x+b)\tan y=0$. additional to (6.4). This
defines another (multi-branched) curve $L$ which has to meet the curve $K$ in
certain points. The function to be considered is

$\ell(x):=-(x+b)\tan k(x).$ (6.7)

The zeros we look for will occur at points where $\ell(x)=k(x)$. It is seen
from Prop. 6.2(a) that this will happen at points where $k(x)$ is near
$(n+1/2)\pi;n=0,1,2,\dots$ and this results in a discrete sequence of zeros of
$h$, each of the form $s=x+ik(x)$.

We next examine the case $x+b\geq 0$ where zeros with non-negative real part
may occur. According to Prop. 6.2(c) this can only happen if $\rho>1$. It will
be necessary to indicate the dependence of functions and derived parameters on
$b=\tilde{\beta}\beta^{{}^{{}_{*}}}=\beta^{{}^{{}_{*}}}\beta\tau$ and
$\rho=\gamma/\beta$. Thus we write $B=B(b,\rho)$, $k(x)=k(x;b,\rho)$ in
accordance with (6.5), (6.6). Also note that $x_{{}_{{}^{M}}}$ also depends on
$b$ and $\rho$.

It is convenient to consider equation (6.2) on the curve $K$ instead of (6.7).
This yields, after some manipulation,

$r(x;b,\rho):=(x+b)\exp\\{x+b\\}+B(b,\rho)\cos k(x;b,\rho)=0.$ (6.8)

To investigate this equation we consider $r(-b;b,\rho)=B(b,\rho)\cos
k(-b;b,\rho)=B(b,\rho)\cos B(b,\rho)$ as can be seen from (6.6). Thus, if
$3\pi/2>B(b,\rho)>\pi/2>1$, $r(-b;b,\rho)<0$. Also at $x=x_{{}_{{}^{M}}}$, we
find that
$r(x_{{}_{{}^{M}}};b,\rho)=(x_{{}_{{}^{M}}}+b)\exp\\{x_{{}_{{}^{M}}}+b\\}+B(b,\rho)>0$
since $k(x_{{}_{{}^{M}}};b,\rho)=0$. We conclude that $r(x;b,\rho)$ has zeros
in the interval $(-b,x_{{}_{{}^{M}}})$. These zeros may still be negative.

We can, however, find a very interesting value of $b$ by noticing that
$k(0;b,\rho)=(\rho^{2}-1)^{-1/2}b$ and

$r(0;b,\rho)=b\exp\\{b\\}[1+\rho\cos([\rho^{2}-1]^{1/2}b)].$ (6.9)

Thus, if we take $b$ as

$b_{{}_{{}^{0}}}=b_{{}_{{}^{0}}}(\rho):=(\rho^{2}-1)^{-1/2}\arccos(-\rho^{-1})=(\rho^{2}-1)^{-1/2}[\pi-\arccos(\rho^{-1})],$

it is seen that $r(0;b_{{}_{{}^{0}}},\rho)=0$. Moreover, since $\rho>1$,
$0<\arccos(\rho^{-1})<\pi/2$ so that
$\pi/2<(\rho^{2}-1)^{1/2}b_{{}_{{}^{0}}}<\pi$. Thus, $\rho
b_{{}_{{}^{0}}}>\pi/2$ and hence
$B(b_{{}_{{}^{0}}},\rho)>(\pi/2)\exp\\{b_{{}_{{}^{0}}}\\}>\pi/2$. For this
particular choice of $b$, $x=0$ therefore is a zero of $r$. Corresponding to
$b_{{}_{{}^{0}}}$ is a critical value of $\tau$:

$\tau_{{}_{{}^{0}}}={b_{{}_{{}^{0}}}}/{\beta\beta^{{}^{{}_{*}}}}=(1+\rho)b_{{}_{{}^{0}}}(\rho)/\beta.$
(6.10)

A positive zero may be contrived by incrementing $b_{{}_{{}^{0}}}$ without
changing $\beta$ and $\gamma$ (i.e., $\rho$ fixed). This amounts to letting
$\tau=\tau_{{}_{{}^{0}}}+\sigma$; $\sigma>0$. Then
$b=b_{{}_{{}^{0}}}+\beta^{{}^{{}_{*}}}\beta\sigma$. From (6.9) we see that if

$\beta\sigma\leq\left[\frac{\rho+1}{\rho-1}\right]^{1/2}[\pi-(\rho^{2}-1)^{1/2}b_{{}_{{}^{0}}}],$

$r(0;b)<0$ which means that a positive zero of (6.8) exists.

The zeros of $h$ are simple. Indeed, if $s^{{}_{{}^{*}}}$ is a zero, then
$s^{{}_{{}^{*}}}+b=-c\exp\\{-s^{{}_{{}^{*}}}\\}$ and
$h^{\prime}(s^{{}_{{}^{*}}})=1-c\exp\\{-s^{{}_{{}^{*}}}\\}=1+b+s^{{}_{{}^{*}}}$
which cannot be zero. Also, when the negative branch of $K$ namely,
$s=x-ik(x)$ is considered, we see from (6.7) that if $s=x+ik(x)$ is a zero on
the positive branch, its complex conjugate $\overline{s}=x-ik(x)$ is a zero on
the negative branch.

Positioning of the zeros of $h$ is illustrated in Fig.1. Increasing $\tau$
could shift $x_{{}_{{}^{0}}}$ to ‘the other side’. We order the roots of $h$
according to their real parts: $x_{{}_{{}^{n}}}>x_{{}_{{}^{n+1}}};n=0,1,\dots$
and note that $x_{{}_{{}^{n}}}\to-\infty$ as $n\to\infty$.

Figure 1: Zeros of $h(s)$

We summarize the findings above:

######  Theorem 6.1.

In terms of the parameters $b=\tilde{\beta}\beta^{{}^{{}_{*}}}$,
$\rho=\gamma/\beta$ and $B(b)=\rho b\exp\\{b\\}$ the following is known about
the zeros of the function $h(s)=s+b[1+\rho\exp\\{-s\\}]$:

  1. (a)

All zeros are simple.

  2. (b)

There is a constant $p>0$ such that all zeros occur to the left of the contour
$C=\\{s=p+iy:y\in\mathbb{R}\\}$.

  3. (c)

If $B(b,\rho)>1$ there is a decreasing unbounded sequence
$\\{x_{{}_{{}^{n}}}<x_{{}_{{}^{M}}}:n=0,1,2,\dots\\}$ such that
$s_{{}_{{}^{n}}}=x_{{}_{{}^{n}}}+ik(x_{{}_{{}^{n}}};b)$ are zeros. The complex
conjugates are also roots.

  4. (d)

If $\rho>1$ and $3\pi/2>B(b,\rho)>\pi/2$, a finite number zeros with non-
negative real part can occur. If $\rho\leq 1$ there are only zeros with
negative real part.

##  7 Laplace inversion

The significance of Thm. 6.1 is seen by considering the Mellin inversion of
$\hat{f}$. From (6.1), $\hat{f}=H(s)/h(s)$ and the inversion yields the
representation

$f(t)=\frac{1}{2\pi i}\int_{C}\frac{H(s)}{h(s)}\exp\\{st\\}\kern 1.00006ptds.$

From the residue theorem the integral on the right equals $2\pi i$ times the
sum of the residues of the integrand at its poles, and the poles are precisely
the zeros of $h$. Now, let
$s_{{}_{{}^{n}}}=x_{{}_{{}^{n}}}+iy_{{}_{{}^{n}}}=x_{{}_{{}^{n}}}+ik(x_{{}_{{}^{n}}})$
be one of the zeros of $h$ obtained from the positive branch of the curve $K$
as discussed so far. To calculate the residue, we consider the Taylor
expansion of $h$ about $s_{{}_{{}^{n}}}$. By taking into account that
$s_{{}_{{}^{n}}}+b=-c\exp\\{-s_{{}_{{}^{n}}}\\}$,
$h^{\prime}(s)=1+b-c\exp\\{-s\\}$ and
$h^{{}^{{}_{(j)}}}(s)=(-1)^{j}c\exp\\{-s\\}$ for $j\geq 2$, we obtain

$h(s)=(s-s_{{}_{{}^{n}}})[1+b+s_{{}_{{}^{n}}}+(b+s_{{}_{{}^{n}}})\sum_{j=2}^{\infty}(-1)^{j-1}\tfrac{(s-s_{{}_{{}^{n}}})^{j-1}}{j!}].$

The pole at $s=s_{{}_{{}^{n}}}$ therefore has the residue

$R(s_{{}_{{}^{n}}})=\tfrac{H(s_{{}_{{}^{n}}})}{1+b+s_{{}_{{}^{n}}}}\exp\\{s_{{}_{{}^{n}}}t\\}.$

However, the negative branch of the curve $K$ also contributes. In fact, if
$s_{{}_{{}^{n}}}=x_{{}_{{}^{n}}}+iy_{{}_{{}^{n}}}$ with
$y_{{}_{{}^{n}}}=k(x_{{}_{{}^{n}}})$, its complex conjugate
$\overline{s}_{{}_{{}^{n}}}$ is also a zero of $h$. The residue at this pole
turns out to be the complex conjugate of $R(s_{{}_{{}^{n}}})$. If we write
$H(s_{{}_{{}^{n}}})/(1+b+s_{{}_{{}^{n}}})=|H(s_{{}_{{}^{n}}})|/|1+b+s_{{}_{{}^{n}}}|\exp\\{i\sigma_{{}_{{}^{n}}}\\}$,
the two residues together contribute to the solution by the term

$f_{{}_{{}^{n}}}(t)=2\frac{|H(s_{{}_{{}^{n}}})|}{|1+b+s_{{}_{{}^{n}}}|}\exp\\{x_{{}_{{}^{n}}}t\\}\cos(y_{{}_{{}^{n}}}t+\sigma_{{}_{{}^{n}}}).$
(7.1)

We immediately note that, since $y_{{}_{{}^{n}}}>0$, there is undulation in
every such term. Also, if $x_{{}_{{}^{n}}}<0$ the term decays to zero
exponentially. This is not the case when $x_{{}_{{}^{n}}}\geq 0$. Thus the
dominant term in the solution will correspond to
$s_{{}_{{}^{0}}}=x_{{}_{{}^{0}}}+ik(x_{{}_{{}^{0}}})$.

We conclude this section by estimating the period $T$ of the principal mode of
undulation, namely that associated with $x_{{}_{{}^{0}}}$. The ‘angular
velocity’ is $y_{{}_{{}^{0}}}=k(x_{{}_{{}^{0}}})$ and the scaled period
$T_{{}_{{}^{\tau}}}=2\pi/y_{{}_{{}^{0}}}$ so that in unscaled time,

$T=T_{{}_{{}^{\tau}}}\tau=\left(\frac{2\pi}{y_{{}_{{}^{0}}}}\right)\tau.$
(7.2)

The value of $x_{{}_{{}^{0}}}$ can be obtained numerically by (carefully)
solving (6.8) for $x$, making sure that the obtained value is the largest.

##  8 Estimates

The aim of this section is to obtain information about long-term behaviour of
the solution of the homogeneous equation $f^{\prime}(t)+bf(t)+cf(t-1)=0$ under
the initial condition $f(t)=\phi(t)$ for $t\in[-1,0];f(0)=f_{{}_{{}^{0}}}$. To
begin with we notice that, at least formally,
$f(t)=\sum_{n=0}^{\infty}f_{{}_{{}^{n}}}(t)$ with the terms $f_{{}_{{}^{n}}}$
given by (7.1). Since $x_{{}_{{}^{n}}}\to-\infty$, there is a smallest $m\geq
0$ such that $x_{{}_{{}^{m}}}<-b$. With is in mind, we define the (possibly)
principal part of the solution as
$f_{{}_{{}^{pr}}}(t):=\sum_{n=0}^{m}f_{{}_{{}^{n}}}(t)$ and the remainder as
$f_{{}_{{}^{rm}}}(t):=\sum_{n=m+1}^{\infty}f_{{}_{{}^{n}}}(t)$ so that
$f(t)=f_{{}_{{}^{pr}}}(t)+f_{{}_{{}^{rm}}}(t)$. We obtain estimates for the
two components under the assumption that $x_{{}_{{}^{0}}}<0$ which means that
all $x_{{}_{{}^{n}}}$ are negative.

First to be considered is the coefficient
$A_{{}_{{}^{n}}}:=2|H(s_{{}_{{}^{n}}})|/|1+b+s_{{}_{{}^{n}}}|$ in (7.1). We
introduce the symbolism $\|\phi\|:=\sup_{t\in[-1,0]}|\phi(t)|$ to obtain from
(6.1)

$\displaystyle|H(s_{{}_{{}^{n}}})|$
$\displaystyle\leq|f_{{}_{{}^{0}}}|+c\exp\\{-x_{{}_{{}^{n}}}\\}\|\phi\|\int_{-1}^{0}\exp\\{-x_{{}_{{}^{n}}}t\\}\kern
1.00006ptdt$
$\displaystyle\qquad\qquad<|f_{{}_{{}^{0}}}|+c|x_{{}_{{}^{n}}}|^{-1}\|\phi\|\exp\\{|x_{{}_{{}^{n}}}|\\}.$
(8.1)

Also,
$|1+b+s_{{}_{{}^{n}}}|^{2}=y^{2}+(x_{{}_{{}^{n}}}+b)^{2}+2(x_{{}_{{}^{n}}}+b)+1>B^{2}\exp\\{-2(x_{{}_{{}^{n}}}+b)\\}+2(x_{{}_{{}^{n}}}+b)$
as can be seen from (6.6). We therefore have

$|1+b+s_{{}_{{}^{n}}}|>c\exp\\{|x_{{}_{{}^{n}}}|\\}\text{ if
}x_{{}_{{}^{n}}}+b\geq 0.$ (8.2)

Cases where $x_{{}_{{}^{n}}}+b<0$ are treated differently. We (temporarily)
set $u_{{}_{{}^{n}}}=-(x_{{}_{{}^{n}}}+b)>0$ to obtain
$|1+b+s_{{}_{{}^{n}}}|^{2}>B^{2}\exp\\{u_{{}_{{}^{n}}}\\}[1-2B^{-2}u_{{}_{{}^{n}}}\exp\\{-2u_{{}_{{}^{n}}}\\}]$.
From the inequality $2u\exp\\{-2u\\}\leq e^{-1}$ we obtain

$|1+b+s_{{}_{{}^{n}}}|>cE\exp\\{|x_{{}_{{}^{n}}}|\\}\text{ if
}x_{{}_{{}^{n}}}+b<0.$ (8.3)

with $E^{2}=1-B^{-2}e^{-1}$. Combination of (8.1) with (8.2) and (8.3) yields

$\displaystyle\tfrac{1}{2}A_{{}_{{}^{n}}}<(c|x_{{}_{{}^{0}}}|)^{-1}\exp\\{-|x_{{}_{{}^{n}}}|\\}|f_{{}_{{}^{0}}}|+\|\phi\|\text{
if }x_{{}_{{}^{n}}}+b\geq 0;$ (8.4)
$\displaystyle\tfrac{1}{2}A_{{}_{{}^{n}}}<E^{-1}[c^{-1}|f_{{}_{{}^{0}}}|\exp\\{-|x_{{}_{{}^{n}}}|\\}+b^{-1}\|\phi\|]\text{
if }x_{{}_{{}^{n}}}+b<0.$ (8.5)

Here use have been made of the inequalities
$|x_{{}_{{}^{n}}}|\geq|x_{{}_{{}^{0}}}|$ and $|x_{{}_{{}^{n}}}|>b$ in the two
different cases.

From (8.4) and (7.1) it is seen that
$|f_{{}_{{}^{pr}}}(t)|<C_{{}_{{}^{pr}}}\exp\\{-|x_{{}_{{}^{0}}}|t\\}$ with
$C_{{}_{{}^{pr}}}$ a positive constant. The infinite series
$f_{{}_{{}^{rm}}}(t)$ needs more attention. For this it is necessary to obtain
information about the behaviour of $x_{{}_{{}^{n}}}$ for $n\geq m+1$. Our
arguments will hinge on the equations (6.2) and (6.4) expressed in the form
$\cos y_{{}_{{}^{n}}}=-B^{-1}(x_{{}_{{}^{n}}}+b)\exp\\{x_{{}_{{}^{n}}}+b\\}$
and
$y_{{}_{{}^{n}}}^{2}=B^{2}\exp\\{-2(x_{{}_{{}^{n}}}+b)\\}-(x_{{}_{{}^{n}}}+b)^{2}$.
Since $u_{{}_{{}^{n}}}:=-(x_{{}_{{}^{n}}}+b)>0$ the equations are

$\displaystyle\cos
y_{{}_{{}^{n}}}=B^{-1}u_{{}_{{}^{n}}}\exp\\{-u_{{}_{{}^{n}}}\\};$ (8.6)
$\displaystyle
y_{{}_{{}^{n}}}^{2}=B^{2}\exp\\{2u_{{}_{{}^{n}}}\\}-u_{{}_{{}^{n}}}^{2}.$
(8.7)

The constant $\epsilon:=[\arccos(B^{-1}e^{-1})]/\pi<1/2$ will provide some
clarity. In §6 it is suggested that $y_{{}_{{}^{n}}}$ should be near
$y^{{}_{{}^{*}}}_{{}_{{}^{n}}}:=(n+1/2)\pi$, and this we shall make more
precise.

From (8.6) we see that $0<\cos
y_{{}_{{}^{n}}}<B^{-1}e^{-1}=\cos(\epsilon\pi)$. Standard trigonometry (even a
good sketch) shows that these inequalities can only be satisfied by
$y_{{}_{{}^{n}}}$ in the intervals
$(y^{{}_{{}^{*}}}_{{}_{{}^{n}}},y^{{}_{{}^{*}}}_{{}_{{}^{n}}}+\epsilon\pi]$
for $n$ odd and
$[y^{{}_{{}^{*}}}_{{}_{{}^{n}}}-\epsilon\pi,y^{{}_{{}^{*}}}_{{}_{{}^{n}}})$
for $n$ even. The endpoints $y_{{}_{{}^{n}}}^{{}_{{}^{*}}}$ are excluded since
$\cos y_{{}_{{}^{n}}}^{{}_{{}^{*}}}=0$.

######  Theorem 8.1.

If $x_{{}_{{}^{n}}}+b<0$ then

$\ln[\exp\\{-(b+1)\\}/(n+1)\pi]<x_{{}_{{}^{n}}}<\ln[c/n\pi].$

###### Proof.

From (8.7) we see that
$y_{{}_{{}^{n}}}^{2}\exp\\{-2u_{{}_{{}^{n}}}\\}=B^{2}[1-(B^{-1}u_{{}_{{}^{n}}}\exp\\{-u_{{}_{{}^{n}}}\\})^{2}]=B^{2}[1-\cos^{2}y_{{}_{{}^{n}}}]=B^{2}\sin^{2}y_{{}_{{}^{n}}}$
(having used (8.6) again). This, in turn gives (after some manipulation)
$x_{{}_{{}^{n}}}=\ln[c|\sin y_{{}_{{}^{n}}}|/y_{{}_{{}^{n}}}]$. Careful
consideration of the cases $n$ even/odd leads to
$B^{-1}e^{-1}=\cos\epsilon\pi\leq|\sin y_{{}_{{}^{n}}}|<1$ and it follows that
$-\ln[\exp\\{b+1\\}y_{{}_{{}^{n}}}\\}\leq
x_{{}_{{}^{n}}}\leq\ln[c/y_{{}_{{}^{n}}}].$ Since $0<\epsilon<1/2$ we have
$y^{{}_{{}^{*}}}_{{}_{{}^{n}}}-\epsilon\pi>n\pi$ and
$y^{{}_{{}^{*}}}_{{}_{{}^{n}}}+\epsilon\pi<(n+1)\pi$. It follows that
$n\pi<y_{{}_{{}^{n}}}<(n+1)\pi$ regardless of the parity of $n$. ∎

The inequality (8.4) may now be employed to estimate $f_{{}_{{}^{n}}}(t)$ in
(7.1) when $n\leq m$. This results in

$|f_{{}_{{}^{n}}}(t)|<2[(c|x_{{}_{{}^{0}}})^{-1}|f_{{}_{{}^{0}}}|\exp\\{-|x_{{}_{{}^{n}}}|\\}+\|\phi\|]\exp\\{-|x_{{}_{{}^{n}}}|t\\}\text{
if }n\leq m.$ (8.8)

For the case $n>m$ we apply (8.5) together with Thm. 8.1 to obtain

$\displaystyle|f_{{}_{{}^{n}}}(t)|$
$\displaystyle<2E^{-1}\left[c^{-1}|f_{{}_{{}^{0}}}|\left(\frac{c}{\pi}\right)^{t+1}\left(\frac{1}{n}\right)^{t+1}\right.$
$\displaystyle\hskip
70.0pt+b^{-1}\|\phi\|\left.\left(\frac{c}{\pi}\right)^{t}\left(\frac{1}{n}\right)^{t}\right]\text{
if }n>m.$ (8.9)

The solution $f$ of the linearized homogeneous problem can now be estimated in
terms of the parameters of the problem:

######  Theorem 8.2.

If $x_{{}_{{}^{0}}}<0$ then:

  1. (a)

There exists a constant $C_{{}_{{}^{pr}}}>0$ such that

$|f_{{}_{{}^{pr}}}(t)|<C_{{}_{{}^{pr}}}\exp\\{-|x_{{}_{{}^{0}}}|t\\}.$

  2. (b)

There exist positive constants $C_{{}_{{}^{rm,1}}},C_{{}_{{}^{rm,2}}}$ such
that for $t>1$

$|f_{{}_{{}^{rm}}}(t)|<C_{{}_{{}^{rm,1}}}\zeta(t+1)\left(\frac{c}{\pi}\right)^{t+1}+C_{{}_{{}^{rm,2}}}\zeta(t)\left(\frac{c}{\pi}\right)^{t}$

with $\zeta$ the Euler-Riemann zeta function.

######  Corollary 8.1.

If $x_{{}_{{}^{0}}}<0$ and $c<\pi$, $f(t)\to 0$ uniformly and exponentially as
$t\to\infty$.

We note that the convergence is determined by $c/\pi$ which factors out of the
partial sums in the series expansion.

##  9 Non-linear behaviour

We now turn to the non-linear equation (5.1), with the usual initial
condition, expressed in the following way:

$\left.\begin{aligned} &f^{\prime}(t)+[bf(t)+cf(t-1)]=F(f(t),f(t-1)),\quad
t>0;\\\ &F(u,v):=-u[bu+cv].\end{aligned}\right\\}$

Our first aim is to obtain a non-linear integral representation. Towards this
we consider the non-homogeneous problem (5.3) once again by
writing$f=f_{{}_{{}^{[0]}}}+f_{{}_{{}^{[v]}}}$ with $f_{{}_{{}^{[0]}}}$ the
solution of the homogeneous problem studied in §6 and $f_{{}_{{}^{[v]}}}$ the
solution of the non-homogeneous problem under the homogeneous initial
condition $f_{{}_{{}^{[v]}}}(t)=0$ for $-1\leq t<0$. For this purpose we
introduce the kernel function $\mathfrak{K}$ as the solution of the problem

$\left.\begin{aligned}
&\mathfrak{K}^{\prime}(t)+b\mathfrak{K}(t)+c\mathfrak{K}(t-1)=0,\text{ for
}t>0;\\\ &\mathfrak{K}(t)=0,\text{ for }t<0;\\\
&\mathfrak{K}(0)=1.\end{aligned}\right\\}$ (9.1)

This differs from the homogeneous problems studied in §§6–8 in the jump
discontinuity at $t=0$ which is not serious. In fact,
$\mathfrak{K}(t)=\exp\\{-bt\\}$ for $0\leq t\leq 1$. Theorems 8.1 and 8.2
apply in this case as well with $f_{{}_{{}^{0}}}=1$ and $\|\phi\|=0$. If
$\lambda:=\ln\\{\pi/c\\}$ this leads to:

$\displaystyle|\mathfrak{K}(t)|$
$\displaystyle<2(c|x_{{}_{{}^{0}}}|)^{-1}\exp\\{-|x_{{}_{{}^{0}}}|(t+1)\\}$
$\displaystyle\hskip 50.0pt+2(Ec)^{-1}\zeta(t+1)\exp\\{-\lambda(t+1)\\}\text{
for }t>1.$

Therefore there are constants $C>0$, $\mu>0$ such that

$|\mathfrak{K}(t)|\leq C\exp\\{-\mu t\\}\text{ for }t\geq 0.$ (9.2)

Repetition of the Laplace transform procedure of §6 for $f_{{}_{{}^{[v]}}}$
yields, in the notation of (6.1),
$\widehat{f_{{}_{{}^{[v]}}}}(s)=\hat{v}(s)/h(s).$

From the convolution theorem we now have
$f_{{}_{{}^{[v]}}}(t)=\int_{0}^{t}\mathfrak{K}(t-t^{\prime})v(t^{\prime})\kern
1.00006ptdt^{\prime}=\int_{0}^{t}\mathfrak{K}(t^{\prime})v(t-t^{\prime})\kern
1.00006ptdt^{\prime},$ and hence,

$f(t)=f_{{}_{{}^{[0]}}}(t)+\int_{0}^{t}\mathfrak{K}(t^{\prime})v(t-t^{\prime})\kern
1.00006ptdt^{\prime}.$ (9.3)

The formal calculations above can be justified by Thm. 5.1 and the (tacit)
assumption that $v$ has a Laplace transform. We continue with the formalism by
letting $v(t)=F(f(t),f(t-1))$.

To see that the integral representation (9.3) is more than formal, we note
that by virtue of Thm.4.2(a) (with $f$, $f_{{}_{{}^{0}}}$ replaced by $f-1$,
$f_{{}_{{}^{0}}}-1$), that the function $t\to v(t)$ is bounded. In fact, if
$|f(t)|\leq M$ then

$|v(t)|\leq\tilde{\beta}M|f(t)|\leq\tilde{\beta}M^{2},$ (9.4)

so that the Laplace transform exists.

It is now possible to consider the asymptotic stability of the nonlinear
equation (5.1) under the initial condition $f(t)=\phi(t)$ for $-1\leq t\leq
0$.

######  Theorem 9.1.

If $\lim_{t\to\infty}v(t)$ exists, $x_{{}_{{}^{0}}}<0$ and $c<\pi$, then
$f(t)-f_{{}_{{}^{[0]}}}\to 0$ as $t\to\infty$.

###### Proof.

From (9.3) we see that

$\displaystyle f_{{}_{{}^{[v]}}}(r+t)-f_{{}_{{}^{[v]}}}(r)$
$\displaystyle=\int_{r}^{r+t}\mathfrak{K}(t^{\prime})v(r+t-t^{\prime})\kern
1.00006ptdt^{\prime}$
$\displaystyle\qquad+\int_{0}^{r}\mathfrak{K}(t^{\prime})[v(r+t-t^{\prime})-v(r-t^{\prime})]\kern
1.00006ptdt^{\prime}.$ (9.5)

To estimate the terms on the right we rely on the inequalities (9.2) and
(9.4).

For the first term on the right of (9.5) one has

$\left|\int_{r}^{r+t}\mathfrak{K}(t^{\prime})v(r+t-t^{\prime})\kern
1.00006ptdt^{\prime}\right|\leq\tilde{\beta}M^{2}C\mu^{-1}\exp\\{-\mu
r\\}[1-\exp\\{-\mu t\\}]$

and this tends to zero as $r\to\infty$.

The second term is treated differently by splitting the integral in two, one
over $[0,r/2]$ and the other over $[r/2,r]$. Since the limit of $v(t)$ exists,
there is for given $\varepsilon>0$, $t_{{}_{{}^{\varepsilon}}}$ such that for
$r,t>t_{{}_{{}^{\varepsilon}}}$, $|v(r+t)-v(r)|<\varepsilon$. Now

$\displaystyle\left|\int_{0}^{r/2}\mathfrak{K}(t^{\prime})[v(r+t-t^{\prime})-v(r-t^{\prime})]\kern
1.00006ptdt^{\prime}\right|$ $\displaystyle\hskip 80.0pt\leq
C\mu^{-1}[1-\exp\\{-\mu r/2\\}]\varepsilon\text{ for
}r,t>2t_{{}_{{}^{\varepsilon}}}.$

The integral over $[r/2,r]$ is treated the same as the first term on the right
of (9.5) and we conclude that $\lim_{t\to\infty}[f(t)-f_{{}_{{}^{[0]}}}(t)]$
exists. The conclusion follows from Cor. 8.1 and Thm. 4.2. ∎

The heuristic argument in §5 for the linear homogeneous equation with solution
$f_{{}_{{}^{[0]}}}$ is justified under the hypotheses of Thm. 9.1. In fact,
since $f_{{}_{{}^{[0]}}}\to 0$ as $t\to\infty$, the theorem shows that under a
fairly weak hypothesis the nonlinear part in the representation (9.3) decays
to zero for large $t$. More explicit conditions for asymptotic stability are
known. In [2, Chap.11] Lyapunov stability of general delay equations is
demonstrated under the assumptions that $x_{{}_{{}^{0}}}<0$ and the initial
state $\phi$ is near the equilibrium level. In [14, Thm. 3] it is shown, by
construction of a Lyapunov function, that asymptotic equilibrium occurs if
$c<b$.

##  10 An example

After the somewhat daunting mathematical sections above, it is appropriate to
give an example. The one presented here is taken from an epidemic which at the
time of this writing was very much on every mind. Local data suggests that the
theta-model (§§2, 3) is necessary. Least squares estimates based on early data
suggest the parameter values $\theta=2.8$, $\beta=0.017$/day and
$f_{{}_{{}^{0}}}=0.046$. The parameters $\tau$ and $\gamma$ have been
manipulated experimentally to obtain results that correspond reasonably to
perceptions. The choices are $\tau=32$ days, $\gamma=0.9$/day. This gives
$B=c\exp\\{b\\}=1.53778<\pi/2$ and $c/\pi=0.475\dots$ All zeros have negative
real part. In fact, a numerical computation based on (6.8) yields
$x_{{}_{{}^{0}}}=-0.043\dots<-b=-0.028\dots$ and $y_{{}_{{}^{0}}}=1.561\dots$

Numerical solution of the initial value problem (3.3), (3.4), with the initial
state taken as constant, namely $\phi\equiv f_{{}_{{}^{0}}}$, (remembering
that $f$ really means $f^{\theta}$) resulted in Fig.2. The period of
undulation, estimated from the linearization, according to (7.2), is $T\approx
2\pi\tau/y_{{}_{{}^{0}}}=128.9$ days. The equilibrium level corresponds to
$F_{{}_{{}^{\infty}}}=[\beta^{{}^{{}_{*}}}]^{1/\theta}=0.238$.

Figure 2: Undulation and decay

The computation reported above shows undulation and suggests decay of
infectivity levels as time goes on. Our analysis of the linear problem in §6
indicates that decay could be to the asymptotic equilibrium state
$f_{{}_{{}^{\infty}}}=1$. We have also shown that by increasing only $\tau$ a
zero with non-negative real part can occur and then, according to the
linearized version, there will be no decay. According to (6.10) the critical
value in the present example is $\tau=\tau_{{}_{{}^{0}}}=34.033\dots$ days.
For this value of $\tau$, $x_{{}_{{}^{0}}}=0$. The question is: would this be
so for the non-linear equation?

First, Thm. 4.2(a) states that infectivity levels cannot run away, but it can
happen in the linearization (according to (7.1)). To come closer to answers it
is instructive to calculate trajectories in the phase portrait ($f^{\prime}$
as a function of $f$). This is shown in Fig.3 for the case discussed above and
with $\tau$ alone increased to 35, slightly above $\tau_{{}_{{}^{0}}}$. The
result seems to confirm that in one case ($\tau=32$) decay is to the
asymptotic point $(1,0)$ and in the other case ($\tau=35$) to a limit cycle
about this point. The linearization indeed leads to some clarification.

Figure 3: Decay to equilibrium (left) and a limit cycle (right)

##  11 The numerical algorithm

A brief note on the method used for the numerical calculations is in order. We
refer to the representation obtained in §4 by the substitution $f=1/g$ and
specifically note the expressions (4.3) and (4.4) used to obtain qualitative
results. For the numerical solution we integrate the latter expression over
the interval $[t,t+\eta]$ to obtain:

$\displaystyle g(t+\eta)$ $\displaystyle=\exp\\{-\delta
i(t)\\}g(t)+\beta^{{}^{{}_{*}}}[1-\exp\\{-\delta i(t)\\}]$
$\displaystyle\qquad+\tilde{\beta}\beta^{{}^{{}_{*}}}\gamma^{{}_{{}^{*}}}\exp\\{-i(t+\eta)\\}\int_{t}^{t+\eta}\exp\\{i(t^{\prime})\\}f(t^{\prime}-1)\kern
1.00006ptdt^{\prime},$

with

$\delta
i(t):=i(t+\eta)-i(t)=\tilde{\beta}[\eta-\gamma^{{}_{{}^{*}}}\int_{t}^{t+\eta}f(t^{\prime}-1)\kern
1.00006ptdt^{\prime}].$

Approximation of the integrals above by the trapezium rule gives

$\displaystyle\int_{t}^{t+\eta}\exp\\{i(t^{\prime})\\}f(t^{\prime}-1)\kern
1.00006ptdt^{\prime}$
$\displaystyle\qquad\qquad\approx\tfrac{1}{2}\eta[\exp\\{i(t+\eta)\\}f(t+\eta-1)+\exp\\{i(t)\\}f(t-1)];$
$\displaystyle\hskip 22.0pt\delta
i(t)\approx\tilde{\beta}\eta\\{1-\tfrac{1}{2}\gamma^{{}_{{}^{*}}}[f(t+\eta-1)+f(t-1)]\\}.$

Combination of everything results in the computational algorithm

$\displaystyle g(t+\eta)$ $\displaystyle\approx g(t)\exp\\{-\delta
i(t)\\}+\beta^{{}^{{}_{*}}}[1-\exp\\{-\delta i(t)\\}]$
$\displaystyle\qquad+\tfrac{1}{2}\eta\tilde{\beta}\beta^{{}^{{}_{*}}}\gamma^{{}_{{}^{*}}}[f(t+\eta-1)+\exp\\{-\delta
i(t)\\}f(t-1)].$

Since $f(t)$ is given for $-1\leq t\leq 0$, this can be computed (coding is
straightforward). The algorithm is grounded in the problem.

##  12 Concluding unscientific remarks

The caption above is borrowed from a similar-sounding title by Johannes
Climacus (Søren Kierkegaard), published in 1846 which, in a way, echoes the
Socratic aphorism that the only wisdom we have is knowing that we do not know
— a paradox that can be (partly) resolved if ‘knowing’ is replaced by
‘understanding’. Understanding, it has been said, expands when horizons of
knowing meet. It is never complete.

A fundamental tenet for the mathematical description of growth is the Verhulst
logistic model which states that growth is determined by what is left to grow
upon. It has the property that growth will increase to devour available
resources. If the initial level is below equilibrium levels will increase
towards equilibrium. On the other hand if the level is initially above
equilibrium it will decrease towards equilibrium. Since
$\beta^{{}^{{}_{*}}}=1-\gamma^{{}_{{}^{*}}}$, the logistic-recovery equation
(3.3) may be re-phrased as

$f^{\prime}(t)=\tilde{\beta}f(t)[1-f(t)]+cf(t)[f(t)-f(t-1)],$

which is the logistic equation perturbed by a recovery term that occurs in Van
der Plank’s equation (1.2) without dormancy. Recovery has the effect of
overshooting the equilibrium state. Logistic growth, acting like a
counterweight, then forces growth to decrease. Once the level is below
equilibrium, growth will turn upwards again. This cyclic process can decay
towards stable equilibrium, but may also become repetitive like the motion of
a pendulum or a planet orbiting the sun.

It is interesting to note that unfettered undulation can only occur when
$\rho=1/R_{{}_{{}^{0}}}>1$ (Thm.6.1). The idea that $R_{{}_{{}^{0}}}<1$ is
prudent, may be questioned. Within the present discussion this can lead to
lowering of the equilibrium level, but at the price of an undulation in which
equilibrium could be a spectre. But then, $R_{{}_{{}^{0}}}$ in the classical
SIR model cannot have the same meaning as $1/\rho$ in the model discussed
here.

The notion of recovery as used here should not be confused with the clinical
use of the word. One might ask: what is recovering, the patient or the
pathogen? The long recovery period (32 days) used in the example makes the
question more incisive; so does the loss of asymptotic equilibrium when
recovery takes longer. One could argue that a sufficiently short recovery
period provides less opportunity for transmission (or evolution) of the
pathogen so that decay to equilibrium would be the result. We should take heed
of the view of Dr. James van der Plank in [16]: observations of infected
subjects merely reflect the state of the pathogen.

Mathematical models serve as a basis for motivated speculation and not much
more. They count among the many metaphors we invent to explain and understand
what is called reality. Computational experiments with such models, not
supported by mathematical insight, are similar to searching for “…two grains
of wheat hid in two bushels of chaff …”. The search for grains of truth may be
long and arduous.

## References

  * [1] R.E. Bellman.  The Stability of Solutions of Linear Differential Equations.  Duke Math. J., 10:643–647, 1943.
  * [2] R.E. Bellman and K.E. Cooke.  Differential-Difference Equations.  Academic Press, 1963.
  * [3] B. Buonomo and M. Cerasuolo.  The Effect of Time Delay in Plant-Pathogen Interactions with Host Demography.  Math Biosc. Eng., 12:473–490, 2015.
  * [4] G. Cacciapaglia, C. Cot, and F. Sannino.  Multiwave pandemic dynamics explained: how to tame the next wave of infectious diseases, March 2021.  In Nature: Scientific Reports, Vol. 11. https://www.nature.com/articles/s41598-021-85875-2.
  * [5] W.J. Cunningham.  A nonlinear difference-differential equation of growth.  Proc. Nat. Acad. Sci., 40:709–713, 1954.
  * [6] M. Della Morte and F. Sannino.  Renormalization group approach to pandemics as a time-dependent SIR model, January 2021.  In Frontiers in Physics, Vol. 8. https://doi.org/10.3389/fphy.2020.591876.
  * [7] W. Feller.  An Introduction to Probability Theory and Its Applications, volume II.  John Wiley & Sons, Inc., 1966.
  * [8] M.E. Gilpin and F.J. Ayala.  Global Models of Growth and Competition.  Proc. Nat. Acad. Sci. USA, 70:3590–3593, 1973.
  * [9] K. Gopalsamy.  Stability and Oscillations in Delay Differential Equations, volume 74 of Mathematics and its Applications.  Springer, Dordrecht, 1992.
  * [10] G.E. Hutchinson.  Circular Causal Systems in Ecology.  Ann. New York Acad. Sci., 50:221–246, 1948.
  * [11] G.S. Jones.  Asymptotic behavior and periodic solutions of a nonlinear differential-difference equation.  Proc. Nat. Acad. Sci., 47:879–882, 1961.
  * [12] W.O. Kermack and A.G. McKendrick.  A Contribution to the Mathematical Theory of Epidemics.  Proc. Roy. Soc. London. Ser. A, 115:700–721, 1927.
  * [13] P.A. Reiser.  Modified SIR Model Yielding a Logistic Solution.  Preprint, 2021.  ArXiv. https://doi.org/10.48550/arXiv.2006.01550.
  * [14] S. Ruan.  Delay Differential Equations in Single Species Dynamics.  In O. Arino, M. Hbid, and E.A. Dads, editors, Delay Differential Equations and Applications, volume 205 of NATO Science Series (II. Mathematics, Physics and Chemistry), pages 477–517. Springer, Dordrecht, 2006.
  * [15] J.E. Van der Plank.  Plant Diseases. Epidemics and Control.  Academic Press, New York, 1963.
  * [16] J.E. Van der Plank.  Dynamics of Plant Disease.  Science, 147:120–124, 1965.
  * [17] P-F Verhulst.  Notice sur la loi que la population suit dans son accroissement.  Correspondance mathématique et physique, 10:113–121, 1838.
  * [18] E.M. Wright.  A nonlinear difference-differential equation.  J. Reine Angew. Math., 194:66–87, 1955.