key: cord-0274418-uvc1pzzm
authors: Berx, Jonas; Indekeu, Joseph O
title: Epidemic processes with vaccination and immunity loss studied with the BLUES function method
date: 2021-06-30
journal: nan
DOI: 10.1016/j.physa.2021.126724
sha: 9b2c4f6a4622ad01dd21d1b160eefc026ae152c7
doc_id: 274418
cord_uid: uvc1pzzm

The Beyond-Linear-Use-of-Equation-Superposition (BLUES) function method is extended to coupled nonlinear ordinary differential equations and applied to the epidemiological SIRS model with vaccination. Accurate analytic approximations are obtained for the time evolution of the susceptible and infected population fractions. The results are compared with those obtained with alternative methods, notably Adomian decomposition, variational iteration and homotopy perturbation. In contrast with these methods, the BLUES iteration converges rapidly, globally, and captures the exact asymptotic behavior for long times. The time of the infection peak is calculated using the BLUES approximants and the results are compared with numerical solutions, which indicate that the method is able to generate useful analytic expressions that coincide with the (numerically) exact ones already for a small number of iterations.

Systems of coupled differential equations (DEs) are omnipresent in the sciences, medicine and even in the humanities. While some of these systems have exact closed-form solutions, the majority do not and hence need to be solved by numerical means. One way to supplement direct numerical simulation with useful information about the (physical) structure behind the coupled DEs is to use approximate (semi-)analytical solutions, in which the "physical" role played by the parameters of the DE is conspicuous.

A stronger need for analytical techniques arises when the parameter space to be scanned is large. For each set of parameters the numerical computation must be repeated, which, for large phase spaces, is computationally expensive. Hence, analytical solutions or approximations are of significant value in this context. Furthermore, numerical techniques can often be unreliable as a result of rounding errors, discretization, etc. For specific problems, one usually requires a specific numerical algorithm that leads to useful results. (Semi-)analytical methods do not rely on the use of discretization, perturbation or linearization and can handle most types of nonlinearities where numerical solvers become very complicated, e.g., integro-differential equations or fractional differential equations.

The best-known examples of methods to calculate such analytical approximants are the Adomian decomposition method (ADM) [1, 2] , the variational iteration method (VIM) [3] and the homotopy perturbation method (HPM) [4] . However, a drawback of these methods is the poor convergence for large values of the independent variable (e.g., for long times).

A well-known example of nonlinear coupled DEs is the susceptible-infected-recoveredsusceptible (SIRS) model for studying the time-evolution of infectious diseases [5] . There have been attempts to generate approximate solutions for this model using the ADM [6] , the VIM and the HPM [7, 8] but these solutions diverge quickly and necessitate an improvement of the methods in order to extend the region of convergence. This paper is structured as follows: in Section II we extend the by now established BLUES function method [9] [10] [11] [12] to a system of nonlinear coupled ordinary DEs. Next, in Section III, we describe the SIRS model with a constant vaccination strategy and reduce the system to the two-dimensional subsystem we subsequently solve. In Section IV we set up the BLUES function method for coupled DEs and apply it, generating approximate solutions to the SIRS model. We stress the importance of the identification of the steady states to construct the associated linear operator. The results are compared with the numerical solution and with results from three other methods, being the ADM, VIM and HPM. We use the analytical approximants to calculate the time at which the infection peak occurs and compare with numerical results. Finally, in Section V we present the conclusions and an outlook.

Here we extend the BLUES iteration originally developed for ordinary DEs [9] [10] [11] and partial DEs [12] to a system of coupled ordinary DEs. The role of the inhomogeneous source (or sink) term in the context of the ordinary DE will now be taken over by a vector of sources (or sinks).

Let us start from an n-dimensional system of inhomogeneous nonlinear coupled ordinary DEs that can be written as a nonlinear operator N t acting on a vector X(t), with source vector χ(t),

We now judiciously decompose the nonlinear operator N t into a linear operator L t , which contains inter alia a first derivative in time, and no higher time derivatives, and a residual operator R, i.e., R ≡ L t −N t , which contains the nonlinear part of N t . We add the subscript t to the linear and nonlinear operators to emphasize their dependence on time. Defining suitable initial conditions for t = 0, i.e.,

completes the description of the system. In the application we consider in this work, R does not depend on t. Thus the action of the linear operator on X results in the following associated linear coupled system

where the subscript denotes derivative w.r.t. time and the same boundary conditions are imposed as in (1) . The elements of the matrix A are constants in most of the applications we have in mind. We now propose to rewrite the system of DEs in an equivalent form by incorporating the initial condition through multiplication of C with a Dirac delta source δ(t) located at t = 0 and by including this term on the right-hand-side of the inhomogeneous system, i.e.,

where we have combined the external source χ and the "initial condition source" Cδ into the combined source ψ. This formulation amounts to resetting the initial condition to zero, so that X(t) = 0 for t ≤ 0 − , followed by a jump in X implied by integrating the DE over the delta source, so that X(0 + ) = C and X(t) evolves in a continuous manner for t > 0.

The solution of this linear system (4) is the following convolution integral [13, 14] 

where G(t) is the Green function matrix for the inhomogeneous linear system. This object can be calculated by finding the matrix exponential exp(At) ≡ 1 + At + ..., i.e.,

with Θ(t) the Heaviside step function (which we take to be unity for positive times including t = 0). This Green function matrix solves the linear system with a delta function unit matrix source, i.e., it is a solution of the matrix equation

Adopting the BLUES function strategy, a solution to the nonlinear system (1), rewritten in the equivalent form

is now proposed in the form of a convolution X(t) = (B * φ)(t), in which the function B(t),

named BLUES function, is taken to be equal to the Green function of the chosen related linear system, i.e., B(t) = G(t) and the new source φ(t) is to be calculated by systematic iteration, using the given (combined) source ψ(t). This procedure starts from the following implicit equation, which makes use of the action of the residual operator,

To find the solution to the nonlinear system (1), equation (9) can be iterated to calculate an approximation for φ in the form of a sequence in powers of the residual R. This leads to the sequence φ (n) (t), with φ (0) (t) = ψ. By subsequently taking the convolution product with B(t), approximate solutions X (n) ψ (t) to (1) can be obtained. Explicitly, the iteration for these approximants reads

where

is the zeroth approximant, which is the convolution product of the linear problem. We now turn to applying the BLUES iteration procedure to the SIRS epidemiological model for the spreading of infectious diseases.

The SIRS model consists of a group of susceptible (S), infected (I) and recovered/immune (R) (human) individuals. The total population N = S + I + R can grow by virtue of a (constant) birth rate π and decay through natural deaths at a rate µ. At birth, the individuals are vaccinated with probability p and consequently acquire immunity, effectively adding up to the group of immune individuals. The remainder of the new births are added to the susceptible group with probability 1 − p. Through contact with infected individuals, a person can become infected, following a mass-action law βSI with force of infection βI. The infected can recover with rate γ, acquiring (temporary) immunity and moving to the group of recovered or immune individuals. Finally, immunity can be lost with rate ξ, whereby recovered individuals move back to the susceptible population. Fig.1 illustrates these different processes in a flow diagram. We assume that all system parameters and populations are positive, and 0 ≤ p ≤ 1.

The interactions between the different populations can be described by the following nonlinear system of DEs, in which the prime denotes derivative w.r.t. time, The time evolution of the total population N = S +I +R can be found by adding (12a)-(12c)

which indicates that the population is not constant. N (t) is a nondecreasing function of t when the birth rate is higher than or equal to the death rate, π ≥ µ. To study the relative importance of the various population fractions, we scale S, I and R by the total population

. This transforms the system (12) into the following system for the fractions,

where s(t) + i(t) + r(t) = 1, ∀t ≥ 0. Note that µ is eliminated by this transformation.

By using the constraint on the population fractions, we can eliminate r(t) and study the "two-dimensional" invariant system

From a stability analysis performed in appendix A, we deduce that the above system has two globally stable fixed points, which are known exactly: a disease-free equilibrium (A1a)

, 0 for which the disease is eradicated and an endemic equilibrium

for which the disease persists and keeps circulating through the population. The final state of the system is characterized by the

which represents the average number of susceptible individuals which are infected by one sick individual during their infectious period, while a vaccination program is in long-time use [15] .

Obviously, the endemic equilibrium ε e can only exist when s * e < 1 and i * e > 0 which means that the vaccination reproduction number must satisfy R V > 1. The disease will be fully eradicated whenever the disease-free equilibrium is the only possible stable fixed point and the endemic equilibrium does not exist. This happens when the vaccination probability p is higher than the critical vaccination threshold p c , which can be inferred from equation

Note that p c may exceed unity, whereas p cannot. Note that one could introduce an active constant vaccination strategy which would amount to moving individuals from S to R with a rate ω. This introduces extra terms into the nonlinear system (12) and results in transforming the term (π + ξ)s(t) into (π + ξ + ω)s(t) in the s−channel of the reduced subsystem (15) . The vaccination reproduction number R V (ω) in this situation would then be equal to

which amounts to rescaling the ω = 0 case, i.e.,

It is now clear that ∀ω > 0, i.e., active vaccination, R V (ω) < R V (0). Hence, active vaccination acts as an additional mechanism to decrease the reproduction number and possibly to suppress the existence of an endemic equilibrium. For the sake of simplicity, we will choose ω = 0 in this work.

To find solutions of the system (15), we can write it as a nonlinear matrix equation, as was demonstrated in section II, i.e.,

with X(t) the vector of solutions,

and with source vector ψ(t) = χ + Cδ(t), with χ the (generally time-dependent) vector of external sources and C the vector of initial conditions, i.e.,

In this work we will choose a constant vaccination strategy at birth, i.e., χ s and χ i are time-independent. Note that we have included the initial conditions in the source ψ by multiplication with a Dirac point source located at t = 0, as was explained in Section II.

Now we judiciously tailor the linear operator that is congruous with the asymptotic equilibrium, by rewriting the nonlinear term in (15) so that already the linear system captures the stable fixed point exactly. This is done by including the deviations of the population fractions from their equilibrium values in the (revised) nonlinear term, as follows,

where s * and i * are the elements of the fixed point vector = (s * , i * ) which represents the equilibrium that is reached. This equilibrium depends uniquely on the value of R V . Note that the refurbished nonlinear term vanishes at the fixed point and represents the product of the fluctuations in susceptible and in infected fractions relative to the equilibrium values.

This approach captures the correct asymptotic behavior for long times provided the linear relaxation times for both s(t) and i(t) exist. This is the case for all R V = 1. For R V = 1 this strategy fails because i(t) is then a marginal variable in the linear system since i (t) = 0 at linear level, which precludes an approach to the fixed point. However, a different choice of linear operator will fix this problem. We will discuss the special (critical) case R V = 1

separately. The novelty of the proposed extension of the BLUES function method to coupled nonlinear systems lies in the above judicious tailoring of the linear system, which includes both equilibria by construction.

With the calibration chosen as in (23) we proceed to identify the linear operator,

where the subscript t on X denotes the time derivative, A is the matrix with elements,

and χ is the vector with elements

The (nonlinear) residual operator R applied to the solution vector X then takes the form

Following the procedure outlined in Section II, we construct an iteration sequence (10) for the solution vector X(t), i.e.,

where B(t) is taken to be the matrix Green function G(t) for the linear problem defined through (24). This G(t) can be found as the inverse of the fundamental matrix of the matrix of coefficients A or equivalently as the matrix exponential of tA multiplied by a step function,

i.e.,

For the disease-free equilibrium (R V < 1) we obtain,

In this (simple) case the Green function matrix is triangular, so that its eigenvalues are conspicuous on the main diagonal. These eigenvalues contain the essential "damping" by virtue of the decaying exponentials, with finite "linear relaxation times" τ s = 1/(π + ξ) and

. Note that τ i diverges for R V ↑ 1. In this limit the relaxation to the disease-free equilibrium becomes "nonlinear". The damping (for R V < 1) ensures that the long-time asymptotics of the approximants, calculated through convolution, are well behaved.

The general Green function matrix, appropriate for both disease-free and endemic equilibria, is more involved and reads,

with

The zeroth approximant (11) is the convolution of the matrix Green function with the source vector ψ(t), i.e.,

and results in the following expressions for the population fraction of susceptible and infected individuals, respectively, 

Upon inspection of this and higher approximants (not reported analytically here) we infer that all BLUES approximants, regardless of the number of iterations (n ≥ 0), are qualitatively correct asymptotically, for all R V = 1, in that they converge exponentially rapidly towards the exact fixed point values for long times, in contrast with the other methods which yield divergences.

We proceed to compare graphically the solution of the SIRS model calculated with the BLUES method with a precise numerical solution and with approximate solutions obtained by the ADM, the VIM, or homotopy perturbation method (HPM). In Table I Table I the equilibrium attained by the system. Case 1: small loss of immunity and high vaccination probability.

As a preliminary remark, we mention that for ξ = 0 (no loss of immunity) the SIRS reduces to a SIR model with vaccination, which was treated earlier in [6, 7] by means of the ADM, HPM and VIM. Here we consider the SIRS model in which the protection offered by vaccination or post-disease immunity is lost with a small probability ξ = 0.1 after some time.

When the vaccination probability p = 0.9 is higher than the critical vaccination threshold p c = 0.5781, the disease will eventually die out and the system will reach the stable diseasefree equilibrium for which i → 0 and s → 0.28. This is shown in Fig. 2 . As we already discussed the BLUES method is accurate and captures the fixed-point values (A1a) of (s * 0 , i * 0 ) in the equilibrium exactly. We remark that the approximants generated by the ADM, VIM and HPM diverge uncontrollably for longer times while the BLUES approximants converge globally for all t ≥ 0 and in every iteration. We also compare the BLUES approximants for different numbers of iteration and notice that they converge rapidly towards the numerical solution. This is shown in Fig. 3 . Case 2: high loss of immunity and high vaccination probability.

As a second example, we consider the case in which immunity is more easily lost (ξ = 0.5) and the population is putting in an effort to vaccinate a larger number of people (p = 0.9).

We can deduce from the critical vaccination probability p c = 1.0406 in Table I that even when all civilians are vaccinated, immunity is lost so quickly that the population always reaches the endemic equilibrium and the disease cannot be eradicated. The result of a comparison between the ADM, VIM, HPM and BLUES methods is shown in Fig. 4 . We also compare the BLUES approximants for different numbers of iteration and observe that they converge rapidly towards the numerical solution. This is shown in Fig. 5 . In this third case, we study the dynamical criticality at R V = 1. The population still reaches the disease-free equilibrium asymptotically, but much more slowly since in the limit R V ↑ 1 the linear relaxation time diverges. The "linear" exponential relaxation is replaced by a "nonlinear" algebraic one, with leading behavior for long times proportional to 1/t. This can be inferred exactly from an analysis of the asymptotic behavior of the system of DEs (23), which at R V = 1 reduces to, with, in this special case, s * = s * 0 = s * e . Inspection of these DEs allows us to establish that the leading asymptotic behavior is a 1/t power-law decay towards the fixed point,

A successful application of the BLUES function method is possible if we acknowledge that we must reconsider the decomposition of the problem into a nonlinear and a linear part. This is necessary in view of the divergence of the "linear relaxation time" τ i = 1/((π +γ)(1−R V )) and the concomitant vanishing of one of the eigenvalues of the linear matrix A given by the triangular form (25), which is appropriate for the disease-free equilibrium. If we stick to this choice there is not enough "damping" in the convolution products that govern the iteration procedure (10) and the approximants diverge, similarly to what routinely happens in the ADM, the VIM and HPM. However, we can recalibrate the linear operator so that both eigenvalues are non-zero and nevertheless the correct (disease-free) fixed point is reached for t → ∞. This is achieved simply by retaining the original nonlinear term in (15) without refurbishing it. The matrix of coefficients A of the linear operator and the residual operator hence become, respectively,

and

This linear system ensures a correct approach to the fixed point for R V = 1 and is a good starting point for the BLUES iteration in this dynamical critical point. We recover global convergence to the numerically exact solution, but at a slower pace than in the noncritical case because the approximants decay exponentially fast for long times, whereas the exact solution features an algebraic decay. The characteristic time of the leading exponentials is constant for all approximants but the amplitude, sign and polynomial prefactors vary as the iteration number n increases. This is how the iteration sequence attempts to approximate an algebraic decay in the limit n → ∞ (see Fig.7 ).

We now calculate the (dimensionless) time t =t at which the peak of the infection occurs from the BLUES approximants and compare with the numerically precise values. At the infection peak, i (t) = 0, and hence, using equation (15b), we deduce that s(t) = (γ + π)/β. implying R V = 1. Note that the BLUES approximant is initially very close to, but later deviates somewhat from the numerical solution. This is due to the difference in the type of asymptotic decay, which is conspicuous in Fig.7 .

model for epidemic spreading. The second is concerned with a technical and methodological refinement, being the development of a matrix BLUES function method for coupled nonlinear ordinary differential equations.

We have shown that the method can be applied to obtain analytic approximants for the SIRS model with vital dynamics, constant vaccination strategy and loss of immunity. We have made a detailed comparison of the iteration procedure with the Adomian decomposition method, the variational iteration method and homotopy perturbation method. It is found that all methods succeed in approximating the (numerically) exact solutions locally and, for the BLUES function method, also globally. The BLUES method generates approximants with a higher accuracy, for the same number of iterations performed, and, moreover, is able to capture the asymptotic behavior of the solutions for long times in both the disease-free and endemic equilibria. In unpublished work [16] we have also studied the SEIRS model, the SIRS model. In general, when using the BLUES function method to generate analytical approximants to systems of coupled DEs, one should tailor the linear system in such a way that it includes all of the existing steady states and at the same time respects the initial condition by means of a source term, in order to achieve globally convergent results.

With the ongoing COVID-19 pandemic, it has become clear that mathematical modeling of epidemic processes is crucial to understanding (and possibly also predicting) the evolution of these viral outbreaks. None of the epidemiological models have exact closed-form solutions, with the exception of the SIS model [17] in which no immunity is possible and individuals return to the susceptible contingent once they recover from the infection. Hence, approximate solutions are needed to assess the impact of the different model parameters without resorting to brute force methods such as numerical simulation. This paper has aimed at showing that the BLUES method is a good candidate for obtaining such approximate analytic solutions. This application of the method to coupled ODEs paves the way for applications to more involved systems such as coupled PDEs [12, 18] , for instance.

of the endemic fixed point (A5) of the shifted system (A4) and hence also of the original system (15) . We can choose the following Lyapunov function [19] V (A7)

Now, from (A4a) and (A4b) it is clear that for the endemic equilibrium the following holds βΣ * e i * e = π(1 − p) − (π + ξ)Σ * e + ξ β (ξ + π) + ξ = (γ + π + ξ)i * e .

(A8)

Substituting this property into (A7) and simplifying gives after some algebra

for all values of Σ, i ≥ 0. This concludes the proof that the endemic equilibrium is globally asymptotically stable. Proving the global stability of the disease-free fixed point is now trivial. One can repeat the previous calculations with the Lyapunov function

which is the same as (A6) with now Σ * 0 instead of Σ * e and i * e → 0.

Solving Frontier Problems of Physics: The Decomposition Method

A review of the decomposition method and some recent results for nonlinear equations

Variational iteration method -Some recent results and new interpretations

Homotopy perturbation technique

A contribution to the mathematical theory of epidemics

Adomian decomposition approach to a SIR epidemic model with constant vaccination strategy

Application of homotopy perturbation and variational iteration methods to SIR epidemic model

Variational iteration and successive approximation methods for a SIR epidemic model with constant vaccination strategy

BLUES function method in computational physics

Analytic iteration procedure for solitons and traveling wavefronts with sources

BLUES iteration applied to nonlinear ordinary differential equations for wave propagation and heat transfer

BLUES function method applied to partial differential equations and analytic approximants for interface growth under shear

Green's matrices for the BVP for the system of ODE of the first order: applications to the beam theories

Ordinary Differential Equations with Constant Coefficients, Translations of mathematical monographs

Reproduction numbers of infectious disease models

Deposition, diffusion and convection: BLUES approximants and some exact results

Three basic epidemiological models

Bacterial range expansions on a growing front: Roughness, fixation, and directed percolation

Lyapunov functions and global stability for SIR, SIRS, and SIS epidemiological models

The system (15) has two fixed points that can be found by a fixed-point analysis which reveals a disease-free equilibrium ε 0 in which the disease has died out, and an endemic equilibrium ε e in which the infected population density reaches a nonzero asymptotic value, i.e.,Note that the disease-free equilibrium ε 0 is independent of the average contact rate β.The endemic equilibrium (A1b) can now be simplified toThese two fixed points are globally stable. This means that for an arbitrary initial condition in the si-plane, the trajectory will converge onto one of these two fixed points. Which fixed point is reached depends on the system parameters π, β, γ, ξ and p.The global asymptotic stability for the endemic equilibrium can be proven as follows.First note that the positive quadrant R 2 + of the si-plane is not an invariant set of the system (15), i.e., when s(t) = 0 then s (t) < 0 for all values i(t) > (π(1 − p) + ξ)/ξ. This can be resolved by shifting (s, i) to (Σ, i), whereHence, the shifted system becomesIt is now easy to see for Σ(t) = 0, now Σ (t) ≥ 0 for all values of i(t). By shifting the system, the endemic equilibrium coordinates change as followsAll global properties of the system remain invariant under the shifting of the coordinates, so we can now try to find a Lyapunov function V (Σ, i) to prove global asymptotic stability