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

The asymptotic stability of the null equilibrium of a linear population model with two physiological structures formulated as a first-order hyperbolic PDE is determined by the spectrum of its infinitesimal generator. We propose an equivalent reformulation of the problem in the space of absolutely continuous functions in the sense of Carath\'eodory, so that the domain of the corresponding infinitesimal generator is defined by trivial boundary conditions. Via bivariate collocation, we discretize the reformulated operator as a finite-dimensional matrix, which can be used to approximate the spectrum of the original infinitesimal generator. Finally, we provide test examples illustrating the converging behavior of the approximated eigenvalues and eigenfunctions, and its dependence on the regularity of the model coefficients.


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., [2,23,28,32,33] and references therein). In particular, we consider the case where individuals are characterized by two different traits [31,43].
In applications to population dynamics, the interest is often focused on the longterm 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 zero equilibrium is determined by the spectrum of the infinitesimal generator (IG) of the semigroup of solution operators [31,43]. 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 [6,7,9,14] and for PDE population models with one structuring variable [4,10,34]. 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 [42]. 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., [16,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 [43,Theorem 3.2]. A similar behavior has been shown for the approximation of R 0 for structured epidemic models [5,8,11], and in the approximation of the solution operators [12,13,15].
For structured models with one single structuring variable, pseudospectral methods have already been proposed to study the stability of the zero equilibrium in [4]. 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, with more structuring variables implementing this technique becomes substantially more involved.
A different approach, which has been successfully employed in the context of nonlinear PDEs with one physiological variable [34] and of renewal equations [35,36], consists in first reformulating the problem at hand via conjugation with an integral operator, and then approximating the resulting transformed operator via pseudospectral (also known as spectral collocation) 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.
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. Finally, we provide some concluding remarks in section 7. 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.
The IG of the semigroup {T(t)} t≥0 is the operator A : D(A) → L 1 (R) defined by where D(A) ⊂ L 1 (R) consists of the functions for which the limit exists. Kang et al. [31,Remark 2.3] prove that the operator A satisfies for a.e. x ∈ [x 0 ,x], a.e. y ∈ [y 0 ,ȳ], and every φ ∈ D(A), and that D(A) satisfies the inclusion while [43,Remark 6.1] claims that equality holds. If φ ∈ 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 zero equilibrium. 10 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 [18,Theorem 9.5] and [24, Theorem VI.1.15]).
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 [31]. However, they can be applied to the case of nontrivial velocities as described in section 6, where we also present an example.

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 [34] in the case of one structuring variable, we argue in terms of the integrated state. In particular, we define an isomorphism between L 1 (R) and a suitable space of functions via integration. We then use this isomorphism and the semigroup {T(t)} t≥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 [39] for further details. For Observe that the double integral in the last term is equal to the iterated integral on the variables a and b in any order, thanks to Fubini's theorem. The space AC(R) of absolutely continuous functions on R in the sense of Carathéodory is a Banach space when equipped with the norm · AC(R) defined as We consider a particular subspace of AC(R), namely Observe that AC 0 (R) is a Banach space, being a closed subspace, and that v( Note that, given a solution u of (2.1)-(2.3), R(x,y) u(t, a, b) da db represents the number of individuals whose structuring variables belong to R(x, y) at time t.
Returning now to (2.1)-(2.3), we define the family of operators {S(t)} t≥0 on AC 0 (R) as S(t) := VT(t)V −1 . Since V and V −1 are linear and bounded, recalling that the composition of a compact and a bounded operator (in either order) is compact, the operators S(t) are in turn linear and bounded and form a family with the same properties as {T(t)} t≥0 , namely they form a strongly continuous and eventually compact semigroup on AC 0 (R). Its IG is B : 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. By [24, Corollary V.3.2(i)], since the semigroups {T(t)} t≥0 and {S(t)} t≥0 are eventually compact, the spectra of both A and B are at most countable and consist only of eigenvalues (of finite algebraic multiplicity). By [15,Proposition 4.1] A and B have the same nonzero eigenvalues (with the same multiplicities) and corresponding eigenvectors. Moreover, observe that A is injective if and only if B is, so 0 is either in both spectra or in none. Therefore, the spectra of A and B coincide. From [31,Proposition 3.1], noting that in that proof

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 Π 0 n,m be the space of bivariate polynomials on R 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 Π 0 n,m as approximations of functions in AC 0 (R). Let Θ X = {x 1 , . . . , x n } be a mesh of n points in (x 0 ,x], with x 0 < x 1 < · · · < x n =x, and let Θ Y = {y 1 , . . . , y m } be a mesh of m points in (y 0 ,ȳ], with y 0 < y 1 < · · · < y m =ȳ. We approximate a function ψ ∈ AC 0 (R) by a vector Ψ ∈ R nm according to where the components of Ψ are ordered according to the lexicographic order of the double indices (i, j). Given The finite-dimensional approximation of the operator B is then B n,m : We can write more explicitly the entries of the matrix B n,m by using the bivariate Lagrange representation of ψ n,m , together with the explicit action of the operator B defined in (3.1). Let { X,i } i=0,...,n and { Y,j } j=0,...,m be the Lagrange bases of polynomials relevant to {x 0 } ∪ Θ X and {y 0 } ∪ Θ Y , i.e., The polynomial ψ n,m can be written as Note that indeed ψ n,m (x, y) = 0 for x = x 0 or y = y 0 . Using (3.1) and (4.1), we get March 8, 2022 Using the linearity of K α and K β , it is easy to characterize the entries of the matrix B n,m . Let D X ∈ R n×n and D Y ∈ R m×m be defined as In other words, D X and D Y are the part of the differentiation matrices associated with {x 0 } ∪ Θ X and {y 0 } ∪ Θ Y , respectively, deleting the first row and the first column. The bivariate differentiation matrices in x and y are where ⊗ denotes the Kronecker product. We can then write where A, B, M ∈ R nm×nm are defined by for k, i = 1, . . . , n and l, j = 1, . . . , m. Note that, if µ is constant, the matrix M is diagonal with diagonal entries equal to µ.
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 ,x] and [y 0 ,ȳ].
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 [42,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 [42,Theorem 8.2]. The two latter properties are often known as spectral accuracy, see [40, chapter 4] and [3, chapter 2]. Furthermore, observe that the relevant differentiation matrices can be computed explicitly [40].
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 [37,Theorem 4.8]. Moreover, it is easy to verify that the Lebesgue constant for the tensor-product Chebyshev extremal nodes in R is the product of the univariate Lebesgue constants in [x 0 ,x] and [y 0 ,ȳ], 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 A, B and M in (4.2)-(4.4) cannot be computed exactly. To approximate the integrals on R we use the Clenshaw-Curtis cubature formula [19], which is based on Chebyshev extremal points and is spectrally accurate [40,41].
To compute the entries of A and B, given a function f defined on [x 0 ,x] and F ∈ R n such that F i = f (x i ), we consider the approximation and similarly for functions defined on [y 0 ,ȳ] (see, e.g., [21]). This approximation can be extended to the double integrals involved in the entries of the matrix M for nonconstant µ. More precisely, given a function ψ ∈ AC 0 (R) and a vector Ψ ∈ R nm such that Ψ k,l = ψ(x k , y l ), each integral x k x 0 y l y 0 ψ(a, b) da db can be approximated by the corresponding (k, l)-th entry of the vector D −1 X D −1 Y Ψ.

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 analytically, 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 Ψ of B n,m provides an approximation ψ n,m of Vφ, where φ is an eigenfunction of A; an approximation of φ is thus given by ∂ x ∂ y ψ n,m . For each example we study the behavior for increasing n = m of the absolute error λ on the known eigenvalue λ and of the absolute error φ in L 1 norm on the known eigenfunction φ, 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 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 possible spurious ones do not affect our analysis.

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 5 Considering Example 1.1, Figure 5.1 shows that the errors reach the machine precision already for n = m = 2. With n = m = 1 the error φ on the eigenfunction is 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 5.1 shows that both λ and φ 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 µ.    coefficients α 1 and β 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. In [31] it is assumed that α and β are Lipschitz continuous in order to prove that the semigroup {T(t)} t≥0 is eventually compact, but actually this hypothesis is used only in the interior of the domains: thus, Example 2.4 still falls inside the scope of [31]. Observe also that the function 1 [0,+∞[ (x − y), even if not continuous, is in D(A). Figure 5.2 suggests that the errors decay with finite order and that these orders increase with the regularity of the eigenfunction. In particular, focusing on Examples 2.1-2.3, we can observe that for both errors a loss of one order of differentiability of the eigenfunction seems to correspond to a loss of about one order of convergence (cfr. the dashed reference lines in Figure 5.2). The convergence of λ seems to be almost two orders faster than that of φ . 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.

Nonsmooth eigenfunctions
As for Example 2.4, it appears that losing continuity of the eigenfunction itself changes the order of convergence differently for λ and φ .

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 [2,23,28,32,33] 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 zero solution hold [43]. 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 [34].
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 A, B and M are the matrices defined in section 4 and G X = G X ⊗ I m , G Y = I n ⊗ G Y , with G X and G Y diagonal matrices defined by [G X ] i,i = g X (x i ), i = 1, . . . , n, and [G Y ] i,i = g Y (y j ), j = 1, . . . , m.
As an example let us consider for which λ = −5 is an eigenvalue of the corresponding IG with eigenfunction φ(x, y) = e x 2 −y 2 . We can observe in Figure 6.1 that the errors computed by our method decay with infinite order even in this case.

Concluding remarks
In this paper we proposed a numerical technique to analyze the stability of the zero solution of linear population models with two structuring variables, of the type considered in [31,43].
Extensive numerical tests illustrate the convergence of the eigenvalues of the finitedimensional 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") [1,17,20]. 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 [1]. 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 [25,26].
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 [7,12,35]. 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., [2,23,30,38]), 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 [1]. 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 [27,36]. It would be interesting to apply similar techniques [3] 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, [4] 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.