A pseudospectral method for investigating the stability of linear population models with two physiological structures

1.   Introduction

Mathematical models are often used to describe the evolution of populations in biology and epidemiology. An important class of models that has attracted increased attention is that of structured population models, in which individuals are characterized by one or more variables that describe the i-state (i.e., the individual state) and determine the individual processes, including for instance birth, growth and death. Examples of physiological structures are age, size, spatial position or time since infection. We here focus on a class of structured population models where the structuring variables are continuous and the models are formulated as first-order hyperbolic partial differential equations (PDEs) (see, e.g., ^[1,2,3,4,5] and references therein). In particular, we consider the case where individuals are characterized by two different traits ^[6,7].

In applications to population dynamics, the interest is often focused on the long-term properties of the systems, for instance the existence of equilibrium states and their stability. For linear models with two structures, it has been proved that the stability of the null equilibrium is determined by the spectrum of the infinitesimal generator (IG) of the semigroup of solution operators ^[8,9]. Since the IG is an operator acting on an infinite-dimensional space of functions, numerical techniques are required to obtain finite-dimensional approximations of the operator and, in turn, of its spectrum.

For the analysis of local stability of equilibria, pseudospectral methods have been widely used both for delay equations ^{[10,11,12,13]} and for PDE population models with one structuring variable ^[14,15,16]. The main advantage of pseudospectral methods is their typical spectral accuracy, by which the order of convergence of the approximation error increases with the regularity of the approximated function. In particular, the convergence is exponential for analytic functions ^[17]. In the case of delay equations, this implies that the spectrum of the IG is approximated with exponential order of convergence, as the corresponding eigenfunctions are exponentials (see, e.g., [18,Proposition 3.4]). In the case of PDEs, the eigenfunctions are still exponential in time, but the order of regularity with respect to the physiological variables depends on the regularity of the model parameters, which therefore affects the order of convergence of the approximation [9,Theorem 3.2]. A similar behavior has been shown for the approximation of $R_0$ for structured epidemic models ^[19,20,21], and the approximation of the solution operators ^[22,23,24].

For structured models with one single structuring variable, pseudospectral methods have already been proposed to study the stability of the null equilibrium in ^[14]. In that paper the IG is approximated by combining pseudospectral differentiation with the inversion of a (linear) algebraic condition characterizing the domain of the operator. However, implementing this technique becomes substantially more involved in the presence of more structuring variables.

A different approach, which has been successfully employed in the context of nonlinear PDEs with one physiological variable ^[16] and of renewal equations ^[25,26], consists in first reformulating the problem at hand via conjugation with an integral operator, and then approximating the resulting transformed operator via pseudospectral techniques. The advantages of this approach are mainly twofold: on the one hand, the transformed operator acts on a space of absolutely continuous functions (rather than the original space $L^1$), hence point evaluation, as well as polynomial interpolation and collocation, are well defined; on the other hand, the domain of the transformed operator is characterized by a trivial condition (specifically, a zero boundary condition), which substantially simplifies the numerical implementation. From a modeling point of view, the integrated state has a clear interpretation as it represents the number of individuals whose i-state is less than a given value.

Goal of this work is to introduce a numerical method for the stability analysis of linear PDE population models with two structuring variables based on this second approach. As far as we know, there is currently no other numerical method available for this problem. We demonstrate the applicability and computational efficiency of the new method with several numerical tests, illustrating the convergence of the approximated eigenvalues to the exact ones and supporting the conjecture of spectral accuracy, and we apply the method to two interesting test problems from epidemiology and ecology.

The paper is organized as follows. In Section 2 we introduce the prototype model and the relevant solution operators and IG. Section 3 describes the reformulation of the IG in terms of the integrated state. The reformulated IG is discretized in Section 4 and the resulting numerical method is applied to some test problems in Section 5. In Section 6 we discuss the extension to models with structuring variables evolving with nontrivial velocities and show an example. In Section 7 we apply the method to an epidemiological model structured by age and time since infection ^[6] and to a population model structured by age and size ^[7]. Finally, we provide some concluding remarks in Section 8. In Appendix A, for completeness, we apply the method to the case of one structuring variable and illustrate it with a test model.

A MATLAB implementation of the method is available at http://cdlab.uniud.it/software.

2.   The prototype model

Let $x_{0}, \bar x, y_{0}, \bar y\in\mathbb{R}$ such that $x_{0} < \bar x$ and $y_{0} < \bar y$ and let ${\Omega} {:=} [x_0, \bar{x}]\times[y_0, \bar{y}]$. We consider the scalar first-order linear hyperbolic PDE

with boundary conditions

where $u(t, x, y)$ is the density of the given population at time $t\geq0$ depending on the two structuring variables $x\in[x_{0}, \bar x]$ and $y\in[y_{0}, \bar y]$. Here, the rate $\mu(x, y)$ represents the per capita mortality rate of individuals in state $(x, y)$, while $\alpha(x, \xi, \sigma)$ and $\beta(y, \xi, \sigma)$ describe the per capita rates at which individuals in state $(\xi, \sigma)$ produce new individuals in state $(x, y_0)$ and $(x_0, y)$, respectively.

Following [8,Assumption 2.1], we assume that the model coefficients $\mu \colon {\Omega} \to \mathbb{R}$, $\alpha \colon [x_0, \bar{x}] \times {\Omega} \to \mathbb{R}$ and $\beta \colon [y_0, \bar{y}] \times {\Omega} \to \mathbb{R}$ are nonnegative functions, with $\mu$ an $L^1_{\text{loc}}$ function^* on $[x_0, \bar{x}) \times [y_0, \bar{y})$, $\alpha$ and $\beta$ $L^1$ functions, $\alpha(x, \cdot, \cdot)$ and $\beta(y, \cdot, \cdot)$ dominated by $L^1$ functions in $x$ and $y$, respectively, and $\mu$ bounded away from $0$ (for a similar but more general approach, see [27,Section 4]). Observe that the operators $K_\alpha$ and $K_\beta$ map $L^1({\Omega})$ to $L^1([x_0, \bar{x}])$ and $L^1([y_0, \bar{y}])$, respectively, and are bounded.

^*In ^[8] $\mu$ is required to be $L^1$ on its domain; here we make a different choice to allow for unbounded mortality (see Section 2.1) ^[28,29]. Accordingly, in order to have a well-posed IG we added a condition on $\mu\phi$ in its domain in (2.6).

Under these assumptions, for every $u_0 \in L^1({\Omega})$ the initial–boundary value problem defined by (2.1)–(2.3) and $u(0, \cdot, \cdot) = u_{0}$ admits a unique solution $u(t, \cdot, \cdot) \in L^1({\Omega})$ for $t \geq 0$. Moreover, the family of solution operators $\{T(t)\}_{t\geq 0}$, defined by $T(t)u_0 = u(t, \cdot, \cdot)$, forms a strongly continuous semigroup of bounded linear operators in the Banach space $L^1({\Omega})$, see [9,Theorem 3.1] or [8,Theorem 2.3].

With the further assumption that $\mu$, $\alpha$ and $\beta$ are Lipschitz continuous on the interior of their domains, Kang et al. prove in [8,Sections 4 and 5.2] that $\{T(t)\}_{t\geq 0}$ is eventually compact and has asynchronous exponential growth.

The IG of the semigroup $\{T(t)\}_{t\geq 0}$ is the operator $A \colon D(A) \to L^1({\Omega})$ defined by

where $D(A) \subset L^1({\Omega})$ consists of the functions for which the limit exists. Kang et al. [8,Remark 2.3] prove that the operator $A$ satisfies

for a.e. $x \in [x_0, \bar{x}]$, a.e. $y \in [y_0, \bar{y}]$, and every $\phi \in D(A)$, and that $D(A)$ satisfies the inclusion

while [9,Remark 6.1] claims that equality holds. If $\phi \in D(A)$ is sufficiently smooth, the action (2.5) of the operator $A$ simplifies and can be expressed as

The spectrum of the IG determines the stability of the null equilibrium.^† More precisely, the latter is asymptotically stable if and only if the spectral abscissa of $A$ is negative and it is unstable if the spectral abscissa is positive (see [30,Theorem 9.5] and [31,Theorem Ⅵ.1.15]).

^†Note that since $A$ is an operator on a real Banach space, in order to define and compute its spectrum the space and the operator need to be complexified. For details see, e.g., [32,Section Ⅲ.7].

Observe that in the model defined by (2.1)–(2.3) the structuring variables evolve with the same velocity as time. In Sections 3 and 4 we present the integral reformulation and the discretization restricting to this special case, as in ^[8]. However, they can be applied to the case of nontrivial velocities as described in Section 6.

2.1.   About unbounded mortality

When dealing with structuring variables with finite spans, it is common to require that the mortality becomes unbounded when approaching the maximum values of the structuring variables. As for the case of models with two physiological structures, by defining

it is reasonable to require for example that

In this case, from a numerical point of view, it can be convenient to argue in terms of $\Pi(x, y){:=} e^{-U(x, y)}$. It is easy in fact to see that $\Pi(x, y)\to 0$ for $x\to \bar x$ or $y\to \bar y$. Moreover, by letting $v(t, x, y){:=}\Pi(x, y)^{-1}u(t, x, y)$, we can observe that if $u$ is a solution of (2.1)–(2.3) with $u(0, x, y) = u_0(x, y)$, then $v(t, x, y)$ solves the PDE

with boundary conditions

and $v(0, x, y) = u_0(x, y)\Pi(x, y)^{-1}$.

3.   Equivalent formulation in a space of absolutely continuous functions

To conveniently handle the boundary conditions in $D(A)$ from a numerical point of view, inspired by the approach of ^[16] in the case of one structuring variable, we argue in terms of the integrated state. In particular, we define an isomorphism between $L^1({\Omega})$ and a suitable space of functions via integration. We then use this isomorphism and the semigroup $\{T(t)\}_{t \geq 0}$ to construct an appropriate semigroup acting on a space of functions with higher regularity. With this goal in mind, we first recall the definition and some properties of absolute continuity in the sense of Carathéodory. We refer the reader to ^[33] for further details.

A function $v$ defined on ${\Omega}$ is absolutely continuous in the sense of Carathéodory if and only if there exist $e_v\in\mathbb{R}$, $f_v\in L^1([x_0, \bar{x}])$, $g_v\in L^1([y_0, \bar{y}])$ and $h_v\in L^1({\Omega})$ such that

Observe that the iterated integral in the last term is equal to the double integral on $[x_0, x] \times [y_0, y]$ and to the iterated integral on the variables $a$ and $b$ in the opposite order, thanks to Fubini's theorem. The space $AC({\Omega})$ of absolutely continuous functions on ${\Omega}$ in the sense of Carathéodory is a Banach space when equipped with the norm $\lVert\cdot\rVert_{AC({\Omega})}$ defined as

We consider a particular subspace of $AC({\Omega})$, namely

Observe that $AC_0({\Omega})$ is a Banach space, being a closed subspace, and that $v(x_0, y) = v(x, y_0) = 0$ for $v\in AC_0({\Omega})$. The operator $V \colon L^1({\Omega}) \to AC_0({\Omega})$ defined by

defines an isomorphism between $L^1({\Omega})$ and $AC_0({\Omega})$, with $V^{-1}\psi = \partial_x \partial_y \psi$ for all $\psi \in AC_0({\Omega})$. Observe that both $V$ and $V^{-1}$ are bounded ($\lVert V\rVert = \lVert V^{-1}\rVert = 1$).

Note that, given a solution $u$ of (2.1)–(2.3), ${\int_{x_0}^x\int_{y_0}^y} u(t, a, b) \mathop{}\!\mathrm{d} b { \mathop{}\!\mathrm{d} a}$ represents the number of individuals whose structuring variables belong to $[x_0, x]\times[y_0, y]$ at time $t$.

Returning now to (2.1)–(2.3), we define the family of operators $\{S(t)\}_{t \geq 0}$ on $AC_0({\Omega})$ as $S(t) {:=} VT(t) V^{-1}$. Since $V$ and $V^{-1}$ are linear and bounded, the operators $S(t)$ are in turn linear and bounded and form a family with the same properties as $\{T(t)\}_{t \geq 0}$, namely they form a strongly continuous semigroup on $AC_0({\Omega})$. Its IG is $B \colon D(B) \to AC_0({\Omega})$, with $B{:=} VAV^{-1}$ and $D(B) = VD(A)$.

With the aim of using $B$ to study the stability properties of (2.1)–(2.3), it is important to understand the relation between the spectra of $A$ and $B$. First of all, from $B = VAV^{-1}$ we can write $B-\lambda I = V(A-\lambda I)V^{-1}$ for any $\lambda\in\mathbb{C}$. This implies that the resolvent sets of $A$ and $B$, and consequently their spectra, coincide. Moreover, $\phi\in L^1(\Omega)$ is an eigenvector of $A$ if and only if $V\phi\in AC_0(\Omega)$ is an eigenvector of $B$ corresponding to the same eigenvalue. Thus, given the invertibility of $V$, $A$ and $B$ also share the same eigenvalues (with the same multiplicity). Eventually, from the statements (ⅰ) and (ⅳ) in [8,Proposition 5.6], as well as the proof of (ⅴ), we can conclude that the spectrum of $A$ only consists of eigenvalues, which implies that the same holds for $B$. From [8,Proposition 3.1], noting that in that proof $\chi = A\phi + \mu \phi$, for $\psi \in D(B)$ we can write

4.   Pseudospectral discretization of the IG

In this Section, we use pseudospectral methods with a tensorial approach to obtain a finite-dimensional approximation of the operator $B$, whose spectrum can be used to determine the stability properties of the system.

For $n$ and $m$ positive integers, let $\Pi_{n, m}^0$ be the space of bivariate polynomials on ${\Omega}$ of degree at most $n$ in the first variable and at most $m$ in the second variable, taking value $0$ at $x = x_0$ and $y = y_0$. These conditions are motivated by the fact that we use polynomials in $\Pi_{n, m}^0$ as approximations of functions in $AC_0({\Omega})$. Let $\Theta_X = \{x_1, \ldots, x_n\}$ be a mesh of $n$ points in $(x_{0}, \bar x]$, with $x_{0} < x_{1} < \dots < x_{n} = \bar x$, and let $\Theta_Y = \{y_1, \ldots, y_m\}$ be a mesh of $m$ points in $(y_{0}, \bar y]$, with $y_{0} < y_{1} < \dots < y_{m} = \bar y$. We approximate a function $\psi \in AC_0({\Omega})$ by a vector $\Psi \in \mathbb{R}^{nm}$ according to

where the components of $\Psi$ are ordered according to the lexicographic order of the double indices $(i, j)$.

Given $\Psi \in \mathbb{R}^{nm}$, let $\psi_{n, m}\in\Pi_{n, m}^0$ be the polynomial interpolating $\Psi$ on $\Theta_X \times \Theta_Y$:

The finite-dimensional approximation of the operator $B$ is then $B_{n, m} \colon \mathbb{R}^{nm} \to \mathbb{R}^{nm}$ defined as

We can write more explicitly the entries of the matrix $B_{n, m}$ by using the bivariate Lagrange representation of $\psi_{n, m}$, together with the explicit action of the operator $B$ defined in (3.1). Let $\{\ell_{X, i}\}_{i = 0, \ldots, n}$ and $\{\ell_{Y, j}\}_{j = 0, \ldots, m}$ be the Lagrange bases of polynomials relevant to $\{x_0\} \cup \Theta_X$ and $\{ y_0\} \cup \Theta_Y$, i.e.,

The polynomial $\psi_{n, m}$ can be written as

Note that indeed $\psi_{n, m}(x, y) = 0$ for $x = x_0$ or $y = y_0$. Using (3.1) and (4.1), we get

Using the linearity of $K_\alpha$ and $K_\beta$, it is easy to characterize the entries of the matrix $B_{n, m}$. Let $D_X\in \mathbb{R}^{n\times n}$ and $D_Y\in \mathbb{R}^{m\times m}$ be defined as

In other words, $D_X$ and $D_Y$ are the part of the differentiation matrices associated with $\{x_0\} \cup \Theta_X$ and $\{ y_0 \} \cup \Theta_Y$, respectively, deleting the first row and the first column. The bivariate differentiation matrices in $x$ and $y$ are

where $\otimes$ denotes the Kronecker product. We can then write

where ${\bf A}, {\bf B}, {\bf M} \in \mathbb{R}^{nm \times nm}$ are defined by

for $k, i = 1, \ldots, n$ and $l, j = 1, \ldots, m$. Note that, if $\mu$ is constant, the matrix ${\bf M}$ is diagonal with diagonal entries equal to $\mu$.

We finally note that, although the matrix $B_{n, m}$ is defined for any set of nodes, the choice of the latter is critical to ensure the convergence of the interpolating polynomials and, in turn, of the elements of the spectrum. In the following numerical experiments, we choose the Chebyshev extremal points in each interval $[x_0, \bar x]$ and $[y_0, \bar y]$.

In the univariate case, these nodes guarantee that the convergence rate of the interpolating polynomial of degree $n$ is $O(n^{-k})$ if the interpolated function is $C^k$ [17,Theorem 7.2], which implies that the order of convergence is infinite if the function is smooth. Moreover, the convergence rate is $O(c^n)$ for some $0 < c < 1$ if the function is analytic [17,Theorem 8.2]. The two latter properties are often known as spectral accuracy, see [34,chapter 4] and [35,chapter 2]. Furthermore, observe that the relevant differentiation matrices can be computed explicitly ^[34].

The classic result on the interpolation error being bounded by means of the best uniform approximation error and the Lebesgue constant holds also in the bivariate case. A multidimensional version of Jackson's theorem on the best uniform approximation error holds as well [36,Theorem 4.8]. Moreover, it is easy to verify that the Lebesgue constant for the tensor-product Chebyshev extremal nodes in ${\Omega}$ is the product of the univariate Lebesgue constants in $[x_0, \bar x]$ and $[y_0, \bar y]$, hence it is $O(\log n \log m)$. The tensor-product Chebyshev interpolation is thus near-optimal also in the bivariate case.

Although a proof of convergence for the method is out of the scope of this paper, we show that the order of convergence observed numerically for the approximated eigenvalues and eigenvectors is consistent with the well-established order of convergence of polynomial interpolation.

For implementation purposes, we observe that in general the integrals defining ${\bf A}$, ${\bf B}$ and ${\bf M}$ in (4.2)–(4.4) cannot be computed exactly. To approximate the integrals on ${\Omega}$ we use the Clenshaw–Curtis cubature formula ^[37], which is based on Chebyshev extremal points and is spectrally accurate ^[34,38].

To compute the entries of ${\bf A}$ and ${\bf B}$, given a function $f$ defined on $[x_0, \bar{x}]$ and $F \in \mathbb{R}^n$ such that $F_i = f(x_i)$, we consider the approximation

and similarly for functions defined on $[y_0, \bar{y}]$ (see, e.g., ^[45]). This approximation can be extended to the double integrals involved in the entries of the matrix ${\bf M}$ for nonconstant $\mu$. More precisely, given a function $\psi \in AC_0({\Omega})$ and a vector $\Psi \in \mathbb{R}^{nm}$ such that $\Psi_{k, l} = \psi(x_k, y_l)$, each integral $\int_{x_0}^{x_k} \int_{y_0}^{y_l} \psi(a, b) \mathop{}\!\mathrm{d} a \mathop{}\!\mathrm{d} b$ can be approximated by the corresponding $(k, l)$-th entry of the vector ${\bf D}_X^{-1}{\bf D}_Y^{-1}\Psi$.

5.   Numerical experiments

In this Section, we present several numerical experiments to investigate how the spectrum of the finite-dimensional operator $B_{n, m}$ approximates the spectrum of $B$, and in turn of $A$, of each problem at hand. For this purpose, we select several parameter sets for which eigenvalues and eigenfunctions of $A$ can be expressed explictly, and we study the convergence of the approximated eigenvalues of $B_{n, m}$ to the analytic ones. As for the eigenfunctions, we stress that, since $B_{n, m}$ represents an approximation of the operator $B$, an eigenvector $\Psi$ of $B_{n, m}$ provides an approximation $\psi_{n, m}$ of $V\phi$, where $\phi$ is an eigenfunction of $A$; an approximation of $\phi$ is thus given by $\partial_x\partial_y \psi_{n, m}$.

For each example we study the behavior for increasing $n = m$ of the absolute error $\epsilon_\lambda$ on the known eigenvalue $\lambda$ and of the absolute error $\epsilon_\phi$ in $L^1$ norm on the known eigenfunction $\phi$, computed via Clenshaw–Curtis cubature.

In all examples we choose

in the boundary conditions (2.2)–(2.3), in order to simplify finding an explicit eigenfunction.

We remark that the parameters are chosen in order to have an analytically known eigenfunction with certain smoothness properties, without regard to any specific biological interpretation.

To compute the spectrum of $B_{n, m}$ we use standard methods (namely MATLAB's $\texttt{eig}$ function). Note that the approximated spectrum may contain spurious eigenvalues (e.g., when $B$ has fewer eigenvalues than the dimension of $B_{n, m}$); however, in our examples we only examine specific eigenvalues, so that the spurious ones do not affect our analysis.

5.1.   Analytic eigenfunctions

We consider a first group of examples presenting an analytic eigenfunction. The choices of the parameters and the resulting eigenvalue and eigenfunction are listed in Table 1. Starting from Example 1.1, where all parameters are constant, we gradually introduce nonconstant coefficients: $\alpha_1$ and $\beta_1$ in Example 1.2, $\gamma$ in Example 1.3 and $\mu$ in Example 1.4.

Considering Example 1.1, Figure 1 shows that the errors reach the machine precision already for $n = m = 2$. With $n = m = 1$ the errors are exactly equal to $0$, which may be explained by the fact that constant functions are interpolated exactly already by polynomials of degree $0$. As $n = m$ increases, the errors increase, possibly due to numerical instability.

Considering now Examples 1.2, 1.3 and 1.4, Figure 1 shows that both $\epsilon_\lambda$ and $\epsilon_\phi$ decay with infinite order. Observe that in Example 1.4 more nodes are needed to reach the error barrier than in Examples 1.2 and 1.3, probably due to the approximation of the integrals involving the nonconstant $\mu$.

5.2.   Nonsmooth eigenfunctions

For the second group of examples, we consider eigenfunctions which are not smooth. The choices of the parameters and the resulting eigenvalue and eigenfunction are listed in Table 2. Observe that the eigenfunctions have the same regularity as the coefficients $\alpha_1$ and $\beta_1$, namely $C^{2}$ for Example 2.1, $C^1$ for Example 2.2 and $C^{0}$ for Example 2.3. We consider also Example 2.4 with discontinuous coefficients and eigenfunction.

Figure 2 suggests that the errors decay with finite order and these orders increase with the regularity of the eigenfunction. In particular, we can observe that for both errors a loss of one order of differentiability, or the loss of continuity, of the eigenfunction seems to correspond to a loss of about one order of convergence (cfr. the dashed reference lines in Figure 2). The convergence of $\epsilon_\lambda$ seems to be almost two orders faster than that of $\epsilon_\phi$. To possibly explain this difference, recall that we are actually collocating the eigenvalue problem for $B$, which means that the eigenvalues are the same as $A$, but the eigenfunctions correspond to integrals of the eigenfunctions of $A$, so the comparison between the eigenfunctions involves differentiating the computed ones.

6.   Structuring variables with nontrivial velocity

We have illustrated the method for systems in which both physiological variables evolve at the same velocity as time. This should not be seen as too restrictive, as systems with more general velocity terms ^[1,2,3,4,5] can in some cases be reduced to (2.1)–(2.3) after a suitable scaling of variables, so that similar theoretical results on the stability of the null solution hold ^[9]. In practice, however, the computation of the change of variables, which in general is defined by the solution of an ODE system, may be expensive, although necessary when the individual parameters (e.g., birth and mortality rates) depend on the original (unscaled) variables. In this case, directly approximating the original problem with nontrivial velocities may be convenient from a computational point of view, as observed in ^[16].

In fact, the transformation via integration can be easily carried out for problems of the form

where the positive functions $g_X(x)$ and $g_Y(y)$ describe the rates of change of $x$ and $y$ in time. In this case, it is straightforward to verify that, given the IG $A$, the operator $B = VAV^{-1}$ admits the representation

which is approximated by a matrix of the form

where ${\bf A}$, ${\bf B}$ and ${\bf M}$ are the matrices defined in Section 4 and ${\bf G}_X = G_X \otimes I_m$, ${\bf G}_Y = I_n \otimes G_Y$, with $G_X$ and $G_Y$ diagonal matrices defined by $[G_{X}]_{i, i} = g_X(x_i)$, $i = 1, \dots, n$, and $[G_{Y}]_{j, j} = g_Y(y_j)$, $j = 1, \dots, m$.

As an example let us consider

for which $\lambda = -5$ is an eigenvalue of the corresponding IG with eigenfunction $\phi(x, y) = e^{x^2-y^2}$. We can observe in Figure 3 that the errors decay with infinite order even in this case.

7.   Applications

7.1.   An epidemiological model structured by age and time since infection

We now consider an epidemic reinfection model proposed in ^[6] with individuals structured by demographic (or chronological) age $a \in [0, a^\dagger]$ and time since infection $\tau \in [0, \tau^\dagger]$. All individuals are assumed to experience a natural mortality $\mu(a)$, and we define the corresponding survival probability

and the average life expectancy of individuals

Susceptible individuals can become infected upon contact with an infectious individual at a rate given by $\beta(\tau)$. Infected individuals recover at a rate $\gamma(\tau)$. After linearization around the disease-free steady state, the linear equation for the infected individuals reads

where $\zeta(t, \tau, a)$ is the density of the infected individuals in the linearization around the disease-free steady state and $w^*(a) = \ell(a)/e_0$ is the age profile of the host population in the demographic steady state.

Following [6,Section 4.2], the reproduction number is given by

with

Observe that $\Gamma(\tau)$ is the probability of an individual not being removed before time since infection $\tau$. Hence we expect the disease-free equilibrium to undergo a bifurcation at $R_0 = 1$ and specifically we expect a real eigenvalue crossing the imaginary axis from left to right.

For the numerical simulations we use parameters as in Table 3, which are simplified parameters but broadly consistent with a short infection like COVID-19. For this choice of parameters, $R_0 = \beta_0 / k$ with

Since discretizing a large interval causes slower convergence of interpolation methods, for the purpose of the numerical experiments we rescale the variables to the interval $[0, 1]$ considering new variables $x {:=} \tau/\tau^\dagger$ and $y {:=} a/a^\dagger$ and $u(t, x, y){:=}\zeta(t, \tau, a)$. The resulting model reads

We further apply to this model the reformulation described in 2.1 to avoid issues with the unboundedness of $\mu$ and $\gamma$.

Figure 4 shows the rightmost part of the spectrum of the IG for values of $R_0$ around the bifurcation at $R_0 = 1$: as anticipated, the rightmost eigenvalue is real and crosses the imaginary axis. The figure shows the eigenvalues of the IG discretized with two different values of $n = m$: the comparison suggests that some of the eigenvalues are spurious, more precisely the ones seemingly not lying on the emerging exponential-like curves (excluding, of course, the spectral abscissa).

7.2.   A population model structured by age and size

To further illustrate the efficacy of the method on realistic models with nontrivial velocity, we consider a model inspired from ecology, with individuals structured by their demographic age $a$ and their size $z$, which grows in time with velocity $g(z)$ ^[7]. We will take parameters inspired by the Daphnia population growth model ^[12,39], but assuming that the algae resource is fixed at a certain value. We consider a slight modification of parameters so that both the survival probability at the maximal age and the growth rate at the maximal size are zero. We consider $a \in [0, a^\dagger]$ and $z \in [z_b, z_m]$ with parameters as in Table 4.

We assume that individuals start reproducing after reaching the reproductive size $z_A$. Their offspring's size is distributed according to a density function $w(z)$ with support in $[z_b, z_A]$. The model for the population density $n(t, a, z)$ reads

For studying the model we can define $\Pi(a) = \exp(-\int_0^a (\mu(\xi)-\mu_0) \mathop{}\!\mathrm{d}\xi)$ and $q(t, a, z){:=}\Pi(a)^{-1}n(t, a, z)$. Then it is easy to see that $q$ satisfies

Moreover, we rescale the age variable to the interval $[0, 1]$ as in Section 7.1.

Figure 5 shows the stability chart for the null equilibrium in the parameters $r_m$ and $\gamma$, while Figure 6 shows the value of the spectral abscissa varying those two parameters separately.

8.   Concluding remarks

In this paper we proposed a numerical technique to analyze the stability of the null solution of linear population models with two structuring variables, of the type considered in ^[8,9].

Extensive numerical tests illustrate the convergence of the eigenvalues of the finite-dimensional approximation to the true eigenvalues of the IG. The numerical tests support the conjecture that the order of convergence of the approximation depends on the regularity of the eigenfunction. A rigorous theoretical proof of the convergence of the approximation is left to future work.

Stability analysis requires not only that the eigenvalues of the IG are approximated accurately, but also that no spurious eigenvalues are to the right of the true spectral abscissa. In fact, in our examples we observe that the approximation of the known eigenvalue is the numerical spectral abscissa, suggesting that the method can be effectively used to study the stability.

Structured population models can also be formulated as renewal equations for the population birth rate (or "recruitment function") ^[40,41,42]. The renewal equation formulation is particularly convenient from the theoretical point of view as one can exploit the principle of linearized stability for nonlinear equations, which can not be proved in general for the PDE formulation ^[40]. Results on the asymptotic behavior of solutions have also been recently proved, under special assumptions, for renewal equations defined on a space of measures, which makes it possible to consider a wider set of solutions compared to PDEs ^[43,44].

It would be interesting to apply pseudospectral methods in the framework of renewal equations (admitting a state space of multivalued functions or even measures), by extending the techniques developed for scalar renewal equations ^[11,22,25]. However, when the structuring variables evolve with nontrivial velocity, the renewal equation formulation requires to explicitly invert the age-structure relation defined implicitly by an ODE system, which suffers from the computational challenges highlighted in Section 6. Hence, as explained therein, directly tackling the PDE formulation may be computationally convenient, as it bypasses the solution of the ODE system.

In this paper we restricted to structuring variables in bounded intervals, as this allows to exploit the highly desirable convergence properties of polynomial interpolation on Chebyshev nodes. However, unbounded domains are common in the modeling literature (e.g., ^[1,2,27,29]), for instance when it is not easy to determine a suitable upper bound for a physiological variable a priori, or when the processes are naturally described by probability distributions with unbounded support (e.g., exponential or Gamma distributions).

In order to numerically treat these problems, truncating the domain can be feasible sometimes, but the accuracy of the approximation would depend on both the size of the truncated domain and the number of nodes in the domain. The latter usually becomes very large because the choice of Chebyshev nodes does not exploit the specific characteristics of the solutions, which usually belong to exponentially weighted spaces ^[40]. For this reason, using exponentially weighted interpolation and Laguerre-type nodes has proved successful and more efficient than domain truncation in the case of delay equations ^[26,46]. It would be interesting to apply similar techniques ^[35] to structured models in the PDE formulation with one or even two structuring variables, although the latter brings in additional complications due to the necessity to rely on multivariate interpolation.

A.   Models with one structuring variable

As recalled in the introduction, ^[14] already provides a pseudospectral method, based on a different approach, to approximate the IG of models with one structuring variable. In this appendix, for completeness, we adapt our approach to this case, providing also an example.

We consider the scalar first-order linear hyperbolic PDE

with boundary condition

As references on single structure models, see ^[28,29,47]; in particular, see [29,Section 1.2] for what concerns this appendix.

If $\mu, \beta \colon [x_0, \bar{x}] \to \mathbb{R}$ are nonnegative functions with $\mu\in L^1_{\text{loc}}([x_0, \bar{x}))$ and $\beta \in L^\infty([x_0, \bar{x}])$, for every $u_0 \in L^1([x_0, \bar{x}])$ the initial–boundary value problem defined by (A.1)–(A.2) and $u(0, \cdot) = u_{0}$ admits a unique solution $u(t, \cdot) \in L^1([x_0, \bar{x}])$ for $t \geq 0$. Moreover, the family of solution operators $\{T(t)\}_{t\geq 0}$, defined by $T(t)u_0 = u(t, \cdot)$, forms a strongly continuous and eventually compact semigroup of bounded linear operators in the Banach space $L^1([x_0, \bar{x}])$. Its IG $A \colon D(A) \to L^1([x_0, \bar{x}])$, defined as in (2.4), can be expressed as

and its domain $D(A)$ can be characterized as

Let us equip $AC([x_0, \bar{x}])$ with the norm $\lVert\cdot\rVert_{AC([x_0, \bar{x}])}$ defined as $\lVert f\rVert_{AC([x_0, \bar{x}])} {:=} \lvert f(x_0)\rvert + \lVert f'\rVert_{L^1([x_0, \bar{x}])}$. Let $AC_0([x_0, \bar{x}])$ be the subspace of $AC([x_0, \bar{x}])$ of functions that are null at $x_0$, which is a Banach space, being a closed subspace. The operator $V \colon L^1([x_0, \bar{x}]) \to AC_0([x_0, \bar{x}])$ defined by $V\phi(x) = \int_{x_0}^x \phi(\sigma) \mathop{}\!\mathrm{d} \sigma$ defines an isomorphism between $L^1([x_0, \bar{x}])$ and $AC_0([x_0, \bar{x}])$, with $V^{-1}\psi = \psi'$ for all $\psi \in AC_0([x_0, \bar{x}])$. Observe that both $V$ and $V^{-1}$ are bounded ($\lVert V\rVert = \lVert V^{-1}\rVert = 1$).

We define the family of operators $\{S(t)\}_{t \geq 0}$ on $AC_0([x_0, \bar{x}])$ as $S(t) {:=} VT(t) V^{-1}$. As in Section 3, we observe that they form a strongly continuous and eventually compact semigroup on $AC_0([x_0, \bar{x}])$ with IG $B \colon D(B) \to AC_0([x_0, \bar{x}])$, with $B{:=} VAV^{-1}$ and $D(B) = VD(A)$. We can also derive the following expression for $B$, given $\psi \in D(B)$:

As in Section 3, we can conclude that $A$ and $B$ have the same spectrum, at most countable and consisting only of eigenvalues (of finite algebraic multiplicity).

As an example, let us choose $[x_0, \bar x] = [0, 2]$, $\mu(x)\equiv 1$ and $\beta(x){:=} e^{-x}$. It can be shown that the only real eigenvalue of the corresponding IG is the unique real solution of the equation

which can be approximated to the machine precision with standard methods (e.g., with MATLAB's $\texttt{fzero}$ we obtain $\lambda = -1.203187869979980$). The relevant eigenfuction is

We can observe in Figure 7 that the errors computed by our method decay with infinite order.

Acknowledgments

We are grateful to Odo Diekmann for suggesting that we refer to the theory of absolute continuity in the sense of Carathéodory.

The authors are members of INdAM Research group GNCS. Davide Liessi and Francesca Scarabel are members of UMI Research group "Modellistica socio-epidemiologica". The work of Simone De Reggi and Davide Liessi was supported by the Italian Ministry of University and Research (MUR) through the PRIN 2020 project (No. 2020JLWP23) "Integrated Mathematical Approaches to Socio-Epidemiological Dynamics" (CUP: E15F21005420006). Francesca Scarabel is supported by the UKRI through the JUNIPER modelling consortium (grant number MR/V038613/1) and gratefully acknowledges the Dame Kathleen Ollerenshaw Trust and the London Mathematical Society for travel grants that facilitated the completion of this project.

Conflict of interest

The authors declare there is no conflict of interest.