Mathematical modeling of population structure in bioreactors seeded with light-controllable microbial stem cells

1.   Introduction

1.1.   Motivation

Mathematical models of disease spreading date back to the beginning of the twentieth century when Kermack and McKendrick published their famous epidemiological SIR model ^[1]. Since its invention, many researchers have relied heavily on these basic assumptions and have established more advances models ^[2,3,4]–only to name a few publications and references therein. Additionally, networks in epidemiology have been recently considered to describe dynamics of disease spreading and spreading patterns ^[5,6,7,8].

Special attention has been currently attracted by structured models which take age or spatial structure into account ^[9,10]. However, transmission rates depend on age structure as well as sex structure in general. For that reason, we develop a simple age- and sex-structured SIR model for short-time prediction because we want to keep modeling as interpretable as possible ^[12]. Therefore, we structure our population by both sexes and same size age groups.

Due to current epidemics like COVID-19 ^[13], we decided to stay with a SIR-typed model because data are suited for this type of models. If we take a closer look at data from Robert-Koch Institute in Germany, the assumption of same size age groups will be acceptable for current data. Theoretically, we have to consider continuous age-structure as presented in ^[10]. After this short motivational introduction, we can state our contributions in this article.

1.2.   Contributions

Our contributions can be summarized as follows.

$1)$ We develop a time-continuous age- and sex-structured SIR model for short-term predictions with time-dependent transmission rates between susceptible and infectious people and time-dependent recovery rates.

$2)$ At first, we show certain properties such as non-negativity and boundedness of solutions.

$3)$ Additionally, we provide a thorough proof of global existence of solutions in time to our proposed system. We need non-negative and boundedness to conclude global existence and global uniqueness of the solution in time from inductive arguments based on Banach's fixed point theorem. This underlines usefulness of fixed point theorems for arguments regarding existence and uniqueness of solutions in different mathematical areas ^[11].

$4)$ Furthermore, we prove monotonicity properties of the global unique solution and investigate analytically that it convergences to a disease-free equilibrium.

$5)$ Afterwards, we introduce a time-discrete problem formulation which heavily relies on an explicit-implicit formulation of the right-hand-side function. As a consequence, our numerical solution scheme becomes unconditionally stable with respect to chosen time increments. We further show that all properties of the time-continuous formulation transfer to the time-discrete case.

$6)$ We finally summarize our numerical solution scheme in pseudo-code and one numerical example stresses our theoretical findings.

1.3.   Structure

Our article is structured as follows. After our motivational introduction of Section 1, we formulate the time-continuous age- and sex-structured SIR model in Section 2. Additionally, we analyze global existence and global uniqueness, non-negativity, boundedness, monotonicity and long-time behavior of the solution of this model. After that, we propose an explicit-implicit numerical solution scheme in Section 3. Here, we show that all properties of our time-continuous model transfer to our time-discrete problem formulation. We present one numerical example to illustrate our theoretical findings in Section 4 and finally, we conclude our article with some remarks on possible future research directions in Section 5.

2.   Time-continuous problem

The aim of this section is the description and analysis of an age- and sex-structured SIR model. For that purpose, we briefly state our model and its assumptions. At first, we prove global existence based on a modified version of Gr{ö}nwall's Lemma. Afterwards, we provide proofs for non-negativity, boundedness, global uniqueness, monotonicity and long-time behavior of our model's solution.

2.1.   Mathematical background material

To especially state global existence and global uniqueness of the solution of our age- and sex-structured SIR model, we need to introduce some theoretical background material regarding nonlinear ordinary differential equations. Let us first recall Lipschitz continuity of a function on Euclidean spaces.

Definition 2.1 ([14,Subsection 3.2]). Let $d_{1}, d_{2} \in \mathbb{N}$. If $S \subset \mathbb{R}^{d_{1}}$, a defined function $\mathbf{F} \colon S \longrightarrow \mathbb{R}^{d_{2}}$ is called Lipschitz continuous on $S$ if there exists a non-negative constant $L \geq 0$ such that

holds for all $\mathbf{x}, \mathbf{y} \in S$. Here, $\lVert \cdot \rVert$ denotes a suitable norm on the corresponding Euclidean space.

Let $U \subset \mathbb{R}^{d_{1}}$ be open, let $\mathbf{F} \colon U \longrightarrow \mathbb{R}^{d_{2}}$. We shall call $\mathbf{F}$ locally Lipschitz continuous if for every point $\mathbf{x_{0}} \in U$ there exists a neighborhood $V$ of $\mathbf{x_{0}}$ such that the restriction of $\mathbf{F}$ to $V$ is Lipschitz continuous on $V$.

We consider an initial-value problem

where $\mathbf{z} \left(t \right) = \left(x_{1} \left(t \right), \ldots, x_{n} \left(t \right) \right)$ denotes our solution vector. Our vectorial function is represented by $\mathbf{G} \left(t, \mathbf{z} \left(t \right) \right) = \left(g_{1} \left(t, \mathbf{z} \left(t \right) \right), \ldots, g_{n} \left(t, \mathbf{z} \left(t \right) \right) \right)$ and $\mathbf{z_{0}} \in \mathbb{R}^{n}$ are our given initial conditions. To conclude global existence, we can apply the following theorem that is a direct consequence of Gr{ö}nwall's lemma.

Theorem 2.2 ([14,Theorem 4.2.1]). If $\mathbf{G} \colon \left[0, \infty \right) \times \mathbb{R}^{n} \longrightarrow \mathbb{R}^{n}$ is locally Lipschitz continuous and if there exist non-negative real constants $B$ and $K$ such that

holds for all $\mathbf{z} \left(t \right) \in \mathbb{R}^{n}$, then the solution of the initial value problem (2.2) exists for all time $t \in \mathbb{R}$ and moreover, it holds

for all $t \in \mathbb{R}$.

Finally, we need Banach's fixed point theorem to derive global uniqueness.

Theorem 2.3 ([15,Theorem V.18]). Let $\left(X, \varrho \right)$ be a complete metric space with the metric mapping $\varrho \colon X \times X \longrightarrow \left[0, \infty \right)$. Let $T \colon X \longrightarrow X$ be a strict contraction, i.e. there exists a constant $K \in \left[0, 1 \right)$ such that $\varrho \left(Tx, Ty \right) \leq K \cdot \varrho \left(x, y \right)$ holds for all $x, y \in X$. Then the map $T$ has a unique fixed point.

2.2.   Time-continuous problem formulation

At first, we define the supremum norm of a continuous function $f \colon \left[0, \infty \right) \longrightarrow \mathbb{R}$. It is given by

An equivalent definition can be given for continuous functions on intervals $\left[a, b \right]$. Let us now state the model's assumptions ^[10,16,17]:

$1)$ The population size $N$ is fixed over time $t$, i.e. $N \left(t \right) = N$ for all $t \in \left[0, \infty \right)$;

$2)$ We divide the population into three homogeneous subetaoups, namely susceptible people (S), infectious people (I) and recovered people (R). We can clearly assign every individual to exactly one subetaoup. Hence, we obtain

for all $t \in \left[0, \infty \right)$;

$3)$ We further distinguish our subetaoups. Let $N_{a} \in \mathbb{N}$ be the number of age groups and let $f$ and $m$ be the subscripts for female and male persons respectively. Let $k \in \left\{ 1, \ldots, N_{a} \right\}$ be arbitrary. We denote the $k$-th female susceptible subetaoup by $S_{f, k}$ and the $k$-th male susceptible subetaoup by $S_{m, k}$. Consequently, it is clear how we denote the infectious and recovered subetaoups;

$4)$ Additionally, no births and deaths occur;

$5)$ The time-varying transmission rates $\beta_{S_{m, j}, I_{s, k}} \colon \left[0, \infty \right) \longrightarrow \left(0, \infty \right)$ are Lipschitz continuous and continuously differentiable for fixed $j \in \left\{ 1, \ldots, N_{a} \right\}$, arbitrary $k \in \left\{ 1, \ldots, N_{a} \right\}$ and arbitrary $s \in \left\{ f, m \right\}$. In addition to that, there exists a positive constant $M_{\beta} > 0$ such that $\lVert \beta_{S_{m, j}, I_{s, k}} \rVert_{\infty} \leq M_{\beta}$ for all $t \geq 0$, arbitrary $s \in \left\{ f, m \right\}$ and arbitrary $j, k \in \left\{ 1, \ldots, N_{a} \right\}$;

$6)$ The time-varying recovery rates $\gamma_{I_{s, k}} \colon \left[0, \infty \right) \longrightarrow \left(0, \infty \right)$ are Lipschitz continuous and continuously differentiable for arbitrary $s \in \left\{ f, m \right\}$ and arbitrary $k \in \left\{ 1, \ldots, N_{a} \right\}$. Additionally, there are positive constants $M_{\gamma} > 0$ and $m_{\gamma} > 0$ such that $\lVert \gamma_{I_{s, k}} \rVert_{\infty} \leq M_{\gamma}$ and $\gamma_{I_{s, k}} \left(t \right) \geq m_{\gamma}$ for all $t \geq 0$, arbitrary $s \in \left\{ f, m \right\}$ and arbitrary $k \in \left\{ 1, \ldots, N_{a} \right\}$.

For abbreviation, we write $g^{\prime} \left(t \right) : = \dfrac{ \text{d} g \left(t \right)}{ \text{d}t}$ for the first derivative of a differentiable function $g$ at time $t$. Our equations of the time-continuous age- and sex-structured SIR model read

with susceptible initial conditions $S_{s, j} \left(0 \right) = S_{1, s, j} > 0$, infectious initial conditions $I_{s, j} \left(0 \right) = I_{1, s, j} \geq 0$ and recovered initial conditions $R_{s, j} \left(0 \right) = R_{1, s, j} \geq 0$ for arbitrary $s \in \left\{ f, m \right\}$ and arbitrary $j \in \left\{ 1, \ldots, N_{a} \right\}$. At least one initial condition of the infectious subetaoups should be positive. Obviously, it holds

such that population size is preserved for all $t \geq 0$.

2.3.   Non-negativity and boundedness

We examine non-negativity and boundedness of (2.6).

Lemma 2.4. We obtain

for arbitrary $s \in \left\{ f, m \right\}$, for all $j \in \left\{ 1, \ldots, N_{a} \right\}$ and for all $t \geq 0$ with respect to (2.6).

Proof. We divide our proof into four parts. Let $s \in \left\{ f, m \right\}$ and $j \in \left\{ 1, \ldots, N_{a} \right\}$ be arbitrary in the following.

1) We consider

since $S_{s, j} \left(t \right)$ is contained in both summands and does not depend on the summation index $k$. Hence, we can put this term outside our considered sum. Division by $S_{s, j} \left(t \right)$ now yields

and since we are able to write $S^{\prime}_{s, j} \left(t \right) = \dfrac{ \text{d} S_{s, j} \left(t \right)}{ \text{d}t}$, we can rewrite this equation by

through separation of variables. By integration on the respective time interval $\left[0, t \right]$, we observe that

holds. We finally obtain

Hence, it holds $S_{s, j} \left(t \right) > 0$ for all $t \geq 0$ by our approach of separation of variables. This procedure is feasible because our initial conditions for susceptible people are positive.

2) We examine

under the initial condition $I_{s, j} \left(0 \right) = I_{1, s, j} \geq 0$ for arbitrary $s \in \left\{ f, m \right\}$ and arbitrary $j \in \left\{ 1, \ldots, N_{a} \right\}$. Let us additionally assume that $I_{s, k} \left(0 \right) = I_{1, s, k} \geq 0$ for arbitrary $s \in \left\{ f, m \right\}$ and arbitrary $k \in \left\{ 1, \ldots, N_{a} \right\}$ with $k \not = j$. At least one initial condition $I_{1, \widetilde{s}, \widetilde{j}}$ should be positive. This implies

for all derivatives of initial conditions for infectious subetaoups where the initial conditions are zero at time $t = 0$ since all $S_{s, j} \left(0 \right) > 0$ by assumption and all $I_{s, k} \left(0 \right) \geq 0$ with at least one positive function $I_{\widetilde{s}, \widetilde{j}} \left(0 \right) > 0$ by assumption. Hence, there exists a time $T_{1} > 0$ such that $I_{s, j} \left(T_{1} \right) > 0$ for all $s \in \left\{ f, m \right\}$ and all $j \in \left\{ 1, \ldots, N_{a} \right\}$. Additionally, it holds $I_{s, j} \left(t \right) \geq 0$ for all $t \in \left[0, T_{1} \right]$ for all $s \in \left\{ f, m \right\}$ and all $j \in \left\{ 1, \ldots, N_{a} \right\}$.

Now, we interpret $T_{1} > 0$ as our new starting point for our argument. We have to distinguish two cases.

Case 1: Let $T_{2} > T_{1}$ and let $I_{s_{1}, j_{1}} \left(T_{2} \right) = 0$ be one function of an infectious subetaoup which is non-negative for all $t \in \left[0, T_{2} \right]$. This is feasible due to continuity of these functions. Let there be at least one function of infectious subetaoups which is positive at $t = T_{2}$. As proven in the previous inequality, this implies $I_{s_{1}, j_{1}}^{\prime} \left(T_{2} \right) > 0$. However, this yields the existence of a positive constant $\delta > 0$ such that $I_{s_{1}, j_{1}} \left(t \right) < 0$ for all $t \in \left(T_{2} - \delta, T_{2} \right)$ by continuity. This contradicts our assumption. Hence, all functions of infectious subetaoups stay non-negative - even positive - in this case. By induction, this even holds on future time subintervals.

Case 2: Let $T_{2} > T_{1}$ and let $I_{s, j} \left(T_{2} \right) = 0$ for all $s \in \left\{ f, m \right\}$ and all $j \in \left\{ 1, \ldots, N_{a} \right\}$. This implies the status of disease-free equilibrium for all future time points.

Hence, (2.6) preserves non-negativity with respect to all infectious subetaoups.

3) By our second property and integration of

on the time interval $\left[0, t \right]$, we obtain

It yields full non-negativity preservation of our non-linear ordinary differential equation system (2.6).

4) Our upper bound is a direct consequence of

for all $t \geq 0$ and our proof is complete.

2.4.   Global existence

We now prove a global existence theorem of (2.6) based on Theorem 2.2.

Theorem 2.5. The non-linear first order ordinary differential equation system (2.6) has at least one global solution, i.e. these possible solutions exist for all $t \geq 0$.

Proof. We define the six vectors

which build our solution vector

Now, we define $\mathbf{G} \colon \left[0, \infty \right) \times \mathbb{R}^{6 \cdot N_{a}} \longrightarrow \mathbb{R}^{N_{a}}$ by (2.6) in a straightforward manner. By applying maximum norms, triangle inequalities, non-negativity and boundedness by Lemma 2.4, we obtain

for all $j \in \left\{ 1, \ldots, N_{a} \right\}$ and this yields

Hence, Theorem 2.2 implies global existence of the system's possible solutions in time.

2.5.   Global uniqueness

Now, we are able to prove global uniqueness of our time-continuous problem formulation (2.6).

Theorem 2.6. The non-linear first order ordinary differential equation system (2.6) has exactly one global unique solution in time.

Proof. 1) At first, we need one inequality for our proof. Let $x_{1}, x_{2}, y_{1}, y_{2} \in \mathbb{R}$ be arbitrary. By the triangle inequality, we obtain

2) Let

be two solutions of our initial value problem (2.6) with same time-varying coefficients and same initial value conditions. Let us consider

for arbitrary $s \in \left\{ f, m \right\}$ and arbitrary $j \in \left\{ 1, \ldots, N_{a} \right\}$. Application of the triangle inequality and assumptions on our time-varying coefficients yields

Since all functions are bounded above by the population size N, we obtain

by application of our inequality from the first step of this proof.

3) Let us now consider

By integration on the time interval $\left[0, \tau \right]$, we obtain

Application of the triangle inequality and the second part of this proof yields

4) Furthermore, it holds

We obtain

5) Combining the previous steps, we conclude

on the time interval $\left[0, \tau \right]$. Choose $\tau : = \dfrac{1}{8 \cdot \left(N_{a} + 1 \right) \cdot \max \left\{ M_{\beta}, M_{\gamma} \right\}}$. This implies

and hence, the solution is unique on the time interval $\left[0, \tau \right]$ by Banach's fixed point theorem. Inductively, all previous steps hold on following time intervals $\left[k \cdot \tau, \left(k + 1 \right) \right]$ with arbitrary $k \in \mathbb{N}$ and initial conditions at time point $t = k \cdot \tau$. Therefore, we conclude that the solution is unique for all $t \geq 0$ which proves our assertion.

2.6.   Monotonicity and long-time behavior

We conclude our analysis of our time-continuous problem formulation (2.6) by an investigation of monotonicity and long-time behavior.

Theorem 2.7. We obtain the following properties for arbitrary $s \in \left\{ f, m \right\}$ and for all $j \in \left\{ 1, \ldots, N_{a} \right\}$:

1) $S_{s, j}$ is monotonically decreasing and there exists a number $S^{\star}_{s, j} \geq 0$ such that $\lim\limits_{t \to \infty} S_{s, j} \left(t \right) = S^{\star}_{s, j}$ holds. Additionally, we obtain $S^{\star}_{s, j} > 0$;

2) $R_{s, j}$ is monotonically increasing and there exists a number $R^{\star}_{s, j} \geq 0$ such that $\lim\limits_{t \to \infty} R_{s, j} \left(t \right) = R^{\star}_{s, j}$;

3) $I_{s, j}$ is Lebesgue-integrable on $\left[0, \infty \right)$ and we get $\lim\limits_{t \to \infty} I_{s, j} \left(t \right) = 0$;

4) Our system (2.6) always converges to a disease-free equilibrium

for all solution functions of (2.6).

Proof. We divide our proof in four parts. Let $s \in \left\{ f, m \right\}$ and $j \in \left\{ 1, \ldots, N_{a} \right\}$ be arbitrary.

1) Since $0 \leq S_{s, j} \left(t \right) \leq N$ and $0 \leq I_{s, j} \left(t \right) \leq N$ hold for all $t \geq 0$ by Lemma 2.4, we obtain $S^{\prime}_{s, j} \left(t \right) \leq 0$ for all $t \geq 0$. By separation of variables, we know that

is valid and this implies

Since $S_{s, j}$ is monotonically decreasing, bounded below by zero and

there exists a positive real number $S^{\star}_{s, j}$ such that we obtain the limit $\lim\limits_{t \to \infty} S_{s, j} \left(t \right) = S^{\star}_{s, j}$.

2) By considering $R^{\prime}_{s, j} \left(t \right) = \gamma_{I_{s, j}} \left(t \right) \cdot I_{s, j} \left(t \right) \geq 0$ from Lemma 2.4, we conclude that $R_{s, j}$ is monotonically increasing. Since $R_{s, j}$ is further bounded above by N according to Lemma 2.4, there exists a positive real number $R^{\star}_{s, j}$ such that $\lim\limits_{t \to \infty} R_{s, j} \left(t \right) = R^{\star}_{s, j}$.

3) We have $R^{\prime}_{s, j} \left(t \right) = \gamma_{I_{s, j}} \left(t \right) \cdot I_{s, j} \left(t \right)$ according to our non-linear differential equation system (2.6). Integration on $\left[0, \infty \right)$ yields

This yields

and hence, $I_{s, j}$ is Lebesgue-integrable on $\left[0, \infty \right)$. This shows $\lim\limits_{t \to \infty} I_{s, j} \left(t \right) = 0$.

4) Remember the notation introduced at the beginning of the proof of Theorem 2.5. By our three aforementioned properties, we obtain the limiting vector

and this vector represents the disease-free equilibrium. Hence, our non-linear differential equation system converges to the disease-free equilibrium. This finishes our proof.

3.   Explicit-implicit numerical solution algorithm

Here, we develop an explicit-implicit time-discrete variant of our time-continuous age- and sex-structured SIR model. We organize this section similar to the previous one. Our constructive goal in this section is to present a numerical solution scheme that captures as many properties of its continuous analogue as possible.

Let us assume that our time interval $\left[0, T \right]$ can be divided by a strictly increasing sequence $\left\{ t_{p} \right\}_{p = 1}^{M}$ for $M \in \mathbb{N}$ with $t_{1} = 0$ and $t_{M} = T$. To distinguish continuous and time-discrete solutions, all time-discrete functions are denoted by $S^{ \text{num}}_{s, j} \left(t_{p} \right)$ for example. We additionally assume that time-continuous and time-discrete time-varying transmission rates and recovery rates coincide for all times.

3.1.   Time-discrete problem formulation

Here, we state our explicit-implicit time-discrete problem formulation

of the time-continuous SIR model (2.6) for all $p \in \left\{ 1, \ldots, M - 1 \right\}$ and for all subscripts of age groups $j \in \left\{ 1, \ldots, N_{a} \right\}$. Our initial conditions read

for arbitrary $s \in \left\{ f, m \right\}$ and all $j \in \left\{ 1, \ldots, N_{a} \right\}$ with at least one initial condition of infectious subetaoups to be positive. For abbreviation, we write in short $\Delta_{p + 1} = \left(t_{p + 1} - t_{p} \right)$ for all $p \in \left\{ 1, \ldots, M - 1 \right\}$ in the following. This explicit-implicit time-discrete problem formulation obviously fulfills

for all $p \in \left\{ 1, \ldots, M - 1 \right\}$.

3.2.   Solvability

Let us proceed with unique solvability of our numerical scheme (3.1).

1) We observe from

that

holds for arbitrary $s \in \left\{ f, m \right\}$, for all $j \in \left\{ 1, \ldots, N_{a} \right\}$ and for all $p \in \left\{ 1, \ldots, M - 1 \right\}$. Here, the sum in the denominator is given by

2) We see from

that

holds for arbitrary $s \in \left\{ f, m \right\}$, for all $j \in \left\{ 1, \ldots, N_{a} \right\}$ and for all $p \in \left\{ 1, \ldots, M - 1 \right\}$.

3) We conclude from

that

holds for arbitrary $s \in \left\{ f, m \right\}$, for all $j \in \left\{ 1, \ldots, N_{a} \right\}$ and for all $p \in \left\{ 1, \ldots, M - 1 \right\}$.

4) Hence, all our computations demonstrate that our numerical solution scheme (3.1) is uniquely solvable. We even infer that, in contrast to typical explicit Euler-time stepping schemes, it is unconditionally stable and we avoid non-linearities as in implicit Euler-time stepping schemes. We summarize our computations and our observations in the following theorem.

Theorem 3.1. Our numerical solution scheme (3.1) is uniquely solvable for all time steps. Additionally, it is also unconditionally stable.

Proof. Follow the above computations in Subsection 3.2.

3.3.   Non-negativity and boundedness

Let us first remark that our initial conditions are non-negative. By induction, it follows that

hold from (3.3) - (3.5) for all $s \in \left\{ f, m \right\}$, all $j \in \left\{ 1, \ldots, N_{a} \right\}$ and all $p \in \left\{ 1, \ldots, M \right\}$. Boundedness is a consequence of (3.2). Thus, we can state the following lemma.

Lemma 3.2. We obtain

for arbitrary $s \in \left\{ f, m \right\}$, for all $j \in \left\{ 1, \ldots, N_{a} \right\}$ and for all $p \in \left\{ 1, \ldots, M \right\}$.

3.4.   Monotonicity and long-time behavior

We continue this section with our theorem on monotonicity and long-time behavior of the solution of our explicit-implicit numerical scheme (3.1).

Theorem 3.3. We have the following properties:

1) The sequence $\left\{ S^{ \text{num}}_{s, j} \left(t_{p} \right) \right\}_{p = 1}^{M}$ is monotonically decreasing and there exists a non-negative real number $S^{\star, \text{num}}$ such that $\lim\limits_{p \to \infty} S^{ \text{num}}_{s, j} \left(t_{p} \right) = S^{\star, \text{num}}$;

2) The sequence $\left\{ R^{ \text{num}}_{s, j} \left(t_{p} \right) \right\}_{p = 1}^{M}$ is monotonically increasing and there exists a non-negative real number $R^{\star, \text{num}}$ such that $\lim\limits_{p \to \infty} R^{ \text{num}}_{s, j} \left(t_{p} \right) = R^{\star, \text{num}}$;

3) The sequence $\left\{ I^{ \text{num}}_{s, j} \left(t_{p} \right) \right\}_{p = 1}^{M}$ fulfills $\lim\limits_{p \to \infty} I^{ \text{num}}_{s, j} \left(t_{p} \right) = I^{\star, \text{num}} = 0$

for arbitrary $s \in \left\{ f, m \right\}$ and for all $j \in \left\{ 1, \ldots, N_{a} \right\}$.

Proof. 1) By Lemma 3.2, we know that the sequence $\left\{ S^{ \text{num}}_{s, j} \left(t_{p} \right) \right\}_{p = 1}^{M}$ is bounded. Again by Lemma 3.2 and (3.3) - (3.5), we get

for arbitrary $s \in \left\{ f, m \right\}$, for all $j \in \left\{ 1, \ldots, N_{a} \right\}$ and for all $p \in \left\{ 1, \ldots, M - 1 \right\}$. Hence, the sequence $\left\{ S^{ \text{num}}_{s, j} \left(t_{p} \right) \right\}_{p = 1}^{M}$ is monotonically decreasing and it thus converges. This implies the existence of a non-negative real number $S^{\star, \text{num}}$ such that $\lim\limits_{p \to \infty} S^{ \text{num}}_{s, j} \left(t_{p} \right) = S^{\star, \text{num}}$ holds.

2) By Lemma 3.2, we know that the sequence $\left\{ R^{ \text{num}}_{s, j} \left(t_{p} \right) \right\}_{p = 1}^{M}$ is bounded. Again by Lemma 3.2 and (3.3) - (3.5), we conclude

for arbitrary $s \in \left\{ f, m \right\}$, for all $j \in \left\{ 1, \ldots, N_{a} \right\}$ and for all $p \in \left\{ 1, \ldots, M - 1 \right\}$. Hence, the sequence $\left\{ R^{ \text{num}}_{s, j} \left(t_{p} \right) \right\}_{p = 1}^{M}$ is monotonically increasing and it thus converges. This yields the existence of a non-negative real number $R^{\star, \text{num}}$ such that $\lim\limits_{p \to \infty} R^{ \text{num}}_{s, j} \left(t_{p} \right) = R^{\star, \text{num}}$ holds.

3) Let us assume the contrary. This implies the existence of a positive real number $I^{\star, \text{num}}$ such that $\lim\limits_{p \to \infty} I^{ \text{num}}_{s, j} \left(t_{p} \right) = I^{\star, \text{num}}$ holds. By (3.4), we then know that all values of the sequence are positive from a certain sequence index. Hence, there exists a positive real number $\tilde{I}^{ \text{num, min}}$ such that $I^{ \text{num}}_{s, j} \left(t_{p} \right) \geq \tilde{I}^{ \text{num, min}}$. Considering

from (3.5), we obtain

and summation by parts yields

from the mentioned time index $L$ as our summation beginning. However, this contradicts our second property. Hence, $\lim\limits_{p \to \infty} I^{ \text{num}}_{s, j} \left(t_{p} \right) = I^{\star, \text{num}} = 0$ holds.

3.5.   Convergence analysis

Here, we want to discuss convergence of our proposed numerical scheme (3.1).

Theorem 3.4. In addition to the assumptions of Subsection 2.2, all solution functions $S_{s, j}, I_{s, j}, R_{s, j} \colon \left[0, \infty \right) \longrightarrow \left[0, N \right]$ are assumed to be continuously differentiable twice with globally bounded first and second derivatives. Additionally, all first derivatives of time-varying transmission rates and time-varying recovery rates are assumed to be globally bounded as well. Let $\Delta_{p} \leq 1$ for all $p \in \mathbb{N}$. If $\max\limits_{p \in \mathbb{N}} \Delta_{p} \to 0$ holds, the discrete solution of the numerical scheme (3.1) converges linearly towards the global unique continuous solution on a considered time interval $\left[0, T \right]$.

Proof. Since this proof become relatively technical, we briefly describe our strategy. At first, local errors between continuous and time-discrete solutions are considered. Afterwards, we need to take into account that errors propagate in time. Finally, we investigate cumulation of these errors which finalizes our proof. We adapt ideas from ^[18] and ^[19]. In general, we follow [19,Satz 74.1] and modify ideas for explicit Eulerian time-stepping schemes because our scheme is a mixture of explicit-implicit parts.

1) For investigation of local errors, we assume that

hold for arbitrary $s \in \left\{ f, m \right\}$ and arbitrary $j \in \left\{ 1, \ldots, N_{a} \right\}$ and we consider the time interval $\left[t_{p}, t_{p + 1} \right]$. Here, we thus only consider one time step and denote solutions by $\widetilde{S_{s, j}^{ \text{num}} \left(t_{p + 1} \right)}$, $\widetilde{I_{s, j}^{ \text{num}} \left(t_{p + 1} \right)}$ and $\widetilde{R_{s, j}^{ \text{num}} \left(t_{p + 1} \right)}$ respectively.

1.1) It first holds

and solving this equation for $\widetilde{S_{s, j}^{ \text{num}} \left(t_{p + 1} \right)}$ yields

We consider $\left| S_{s, j} \left(t_{p + 1} \right) - \widetilde{S_{s, j}^{ \text{num}} \left(t_{p + 1} \right)} \right|$. It holds

Zero addition and application of the triangle inequality implies

We define the two terms

and

For $I_{a}$, we obtain

Application of the mean value theorem of calculus yields the existence of $\xi_{a} \in \left(t_{p}, t_{p + 1} \right)$ such that

holds. This implies

For $I_{b}$, we obtain

Application of the triangle inequality and rearranging yields

Since

is valid, we obtain

By the mean value theorem of calculus, there exists $\xi_{b} \in \left(t_{p}, t_{p + 1} \right)$ such that

holds. This implies

Here, ${ \beta ^{\prime}}$ denotes the vector of all derivatives of time-varying transmission rates. We conclude

and summarizing our results, this implies

1.2) From

we obtain

We consider $\left| I_{s, j} \left(t_{p + 1} \right) - \widetilde{I_{s, j}^{ \text{num}} \left(t_{p + 1} \right)} \right|$. It holds

Rearranging of these terms and application of the triangle inequality yields

We define the following three terms

and

$I_{c}$ can be rewritten as

By the mean value theorem of calculus, there exists $\xi_{c} \in \left(t_{p}, t_{p + 1} \right)$ such that

holds. This implies

For $I_{d}$, we obtain

Application of the triangle inequality implies

By the mean value theorem of calculus, there is $\xi_{d} \in \left(t_{p}, t_{p + 1} \right)$ such that

holds. Hence, we conclude

Here, ${ \gamma ^{\prime}}$ denotes the vector containing all derivatives of time-varying recovery rates. We consider

By zero addition, we obtain

Application of the triangle inequality yields

We define the following three terms

and

Considering $I_{e, 1}$, there exists $\xi_{e, 1} \in \left(t_{p}, t_{p + 1} \right)$ such that

holds due to the mean value theorem of calculus. Hence, we obtain

By the mean value theorem of calculus, there exists $\xi_{e, 2} \in \left(t_{p}, t_{p + 1} \right)$ such that

is valid. Application of the triangle inequality yields

Now, we consider

By application of the triangle inequality, we obtain

By inequality (3.6) from Step 1.1), we know that

holds. This implies

Combining our results, we obtain

We define

We conclude

1.3) It holds

We consider $\left| R_{s, j} \left(t_{p + 1} \right) - \widetilde{R_{s, j}^{ \text{num}} \left(t_{p + 1} \right)} \right|$ and obtain

Application of zero addition and the triangle inequality yields

By the mean value theorem of calculus, there are $\xi_{f, 1}, \xi_{f, 2}, \xi_{f, 3}, \xi_{f, 4} \in \left(t_{p}, t_{p + 1} \right)$ such that

hold. This implies

and

By inequality (3.7) from Step 1.2), we know that

is valid. We infer that

holds. Summarizing our results, we obtain

and

in a short manner.

1.4) Conclusively, we obtain

from the inequalities (3.6), (3.7) and (3.8).

2) In reality, the points $\left(t_{p}, S_{s, j}^{ \text{num}} \left(t_{p} \right) \right)$, $\left(t_{p}, I_{s, j}^{ \text{num}} \left(t_{p} \right) \right)$ and $\left(t_{p}, R_{s, j}^{ \text{num}} \left(t_{p} \right) \right)$ do not lie on the continuous solution graph. For that reason, we must investigate how procedural errors $S_{s, j}^{ \text{num}} \left(t_{p} \right) - S_{s, j} \left(t_{p} \right)$, $I_{s, j}^{ \text{num}} \left(t_{p} \right) - I_{s, j} \left(t_{p} \right)$ and $R_{s, j}^{ \text{num}} \left(t_{p} \right) - R _{s, j} \left(t_{p} \right)$ propagate to the $\left(p + 1 \right)$-th time step. These estimates are going to be carried out in the following steps 2) and 3) of this proof.

2.1) At first, we must consider $\left| S_{s, j}^{ \text{num}} \left(t_{p + 1} \right) - \widetilde{S_{s, j}^{ \text{num}} \left(t_{p + 1} \right)} \right|$. Remember that $\widetilde{S_{s, j}^{ \text{num}} \left(t_{p} \right)} = S_{s, j} \left(t_{p} \right)$. Note that

and

are valid. Hence, we obtain

Application of the triangle inequality and zero addition yields

Summarizing this result, we obtain

2.2) Now, we consider $\left| I_{s, j}^{ \text{num}} \left(t_{p + 1} \right) - \widetilde{I_{s, j}^{ \text{num}} \left(t_{p + 1} \right)} \right|$. We first observe that

and

are valid from step 1.2). Application of the triangle inequality and zero addition yields

Using (3.10), we obtain

and the result reads

2.3) We consider $\left| R_{s, j}^{ \text{num}} \left(t_{p + 1} \right) - \widetilde{R_{s, j}^{ \text{num}} \left(t_{p + 1} \right)} \right|$. From step 1.3), we know that

and

hold. By application of the triangle inequality, this implies

Using inequality (3.11), we obtain

We conclude that

holds.

2.4) Now, we want to combine our results. Since $s \in \left\{ f, m \right\}$ and $j \in \left\{ 1, \ldots, N_{a} \right\}$ are arbitrary indices, we infer by inequalities (3.10), (3.11) and (3.12) that

holds because $\Delta_{p + 1} \leq 1$ by assumption. This yields

3) Finally, we can finish our proof of convergence. For abbreviation, we write

where $\mathbf{z} \in \mathbb{R}^{6 \cdot N_{a}}$ is defined as in the proof of Theorem 2.5. Our proof is heavily based on the inequality

for all $x \geq 0$. Note that $t_{1} = 0 < t_{2} < \ldots < t_{M - 1} < t_{M} = T$.

3.1) At first, we want to inductively prove that

holds for all $p \in \left\{ 0, \ldots, M - 1 \right\}$. Let $p = 0$ first. The inequality (3.14) is fulfilled. Let $p = 1$ to understand the concept. By application of the triangle inequality and inequalities (3.9) and (3.13), we see that

is valid. We now assume that

holds. We now want to show that (3.14) follows. We see that

holds. This proves (3.14) by induction.

3.2) Concluding our proof, we consider

from (3.14). We define $\Delta : = \max\limits_{r \in \left\{ 2, \ldots, M \right\}} \Delta_{r}$. We infer that

holds. If the initial conditions of our continuous and our time-discrete problem formulation coincide and $\Delta \to 0$, the time-discrete solution convergences linearly towards the continuous solution. This proves our assertion.

3.6.   Numerical solution algorithm

We briefly summarize our numerical solution algorithm for the time-discrete explicit-implicit numerical scheme (3.1) in Table 1. This summary is intended to give a brief overview of aspects which need to be considered during implementation. Especially, we state all inputs which are important for our time-discrete numerical scheme.

4.   Numerical example

In this section, we illustrate our theoretical findings by one synthetic data example. At first, we sum up all important information to set calculations up. Finally, we show the results and discuss these findings with respect to our theoretical results.

4.1.   Setting

Let us provide our setting. In Table 2, we summarize the corresponding indices of population subetaoups. The total population is divided into six subetaoups. Now, we report the (time-varying) transmission rates $\beta_{S_{s, j}, I_{q, k}} \colon \left[0, \infty \right) \longrightarrow \left[0, \infty \right)$ and (time-varying) recovery rates $\gamma_{I_{s, j}} \colon \left[0, \infty \right) \longrightarrow \left[0, \infty \right)$ for arbitrary $s, q \in \left\{ f, m \right\}$ and arbitrary $j, k \in \left\{ 1, \ldots, N_{a} \right\}$. These data can be found in Tables 3 and 4. This is an imaginary disease which spreads mainly in the adult population. All initial conditions of populations subetaoups are described in Table 5. The final time is set $T = 180$ with an equidistant time sequence

and this implies $M = 181$. The total population size reads $N = 100000$ due to Table 5. Hence, all data are available for our numerical simulation.

4.2.   Results

Here, we present the results of our setting described before. In Figure 1, the temporal development of all susceptible population subetaoups is depicted. It can be clearly seen that the resulting graphs are decreasing in time. In Figure 2, all graphs of the temporal development with regard to all infectious subetaoups are portrayed. Figure 3 illustrates the temporal development of all recovered population subetaoups. As expected, these curves are increasing in time. Finally, conservation of the total population size for our implicit-explicit numerical solution scheme is shown in Figure 4.

5.   Conclusion and outlook

We introduced an age- and sex-structured SIR model for short-term predictions in Section 2. We established global existence, global uniqueness, non-negativity and boundedness of the solution. Additionally, we showed some monotonicity properties and proved convergence to a disease-free equilibrium in the continuous setting. Afterwards, we proposed an explicit-implicit numerical solution scheme in Section 3. We were able to demonstrate that all aforementioned properties transfer to this time-discrete setting of the age- and sex-structured SIR model for short-term predictions. We also concluded that this scheme is linearly convergent towards the continuous solution. For short-term predictions, effects of demography and transmission between age groups can be simplified or neglected in this case.

To continue this work and extend it to long-term predictions that definitely play an important role, it might be fruitful to additionally take birth rates and death rates into account. The works ^[20,21] can serve as examples for extensions of our work. Incubation times also lead to delays from transfer between different compartments. Hence, introduction of delays in our system might be another possible future research direction. Examples can be seen in ^[22]. Furthermore, spatial inhomogeneities should also be considered because spreading of diseases differ in regions depending on social status for example ^[23,24], which leads to ODE-PDE coupled systems. Application of higher-order methods might be considerable as well ^[25,26].

Finally, we stress the fact that the inverse problem in dynamics of biological systems needs further investigation ^{[27,28,29,30]}.

Authors' contributions

Both authors conceived and designed the research. Benjamin Wacker analyzed the time-continuous problem formulation. Benjamin Wacker analyzed the time-discrete problem formulation. Benjamin Wacker implemented the explicit-implicit numerical solution scheme. Both authors discussed the numerical example. Both authors drafted and edited this manuscript.

Conflict of interest

Both authors declare that they have no conflict of interest.