Perturbation analysis in finite LD‐QBD processes and applications to epidemic models

Summary In this paper, our interest is in the perturbation analysis of level‐dependent quasi‐birth‐and‐death (LD‐QBD) processes, which constitute a wide class of structured Markov chains. An LD‐QBD process has the special feature that its space of states can be structured by levels (groups of states), so that a tridiagonal‐by‐blocks structure is obtained for its infinitesimal generator. For these processes, a number of algorithmic procedures exist in the literature in order to compute several performance measures while exploiting the underlying matrix structure; among others, these measures are related to first‐passage times to a certain level L(0) and hitting probabilities at this level, the maximum level visited by the process before reaching states of level L(0), and the stationary distribution. For the case of a finite number of states, our aim here is to develop analogous algorithms to the ones analyzing these measures, for their perturbation analysis. This approach uses matrix calculus and exploits the specific structure of the infinitesimal generator, which allows us to obtain additional information during the perturbation analysis of the LD‐QBD process by dealing with specific matrices carrying probabilistic insights of the dynamics of the process. We illustrate the approach by means of applying multitype versions of the susceptible‐infective (SI) and susceptible‐infective‐susceptible (SIS) epidemic models to the spread of antibiotic‐sensitive and antibiotic‐resistant bacterial strains in a hospital ward.


A MATHEMATICAL MODEL OF BACTERIAL TRANSMISSION
We link the SI 1 , I 2 and SI 1 , I 2 S epidemic models to the deterministic model in [2, Figure A] for the spread of two bacterial strains in a hospital ward. Lipsitch et al. [2] consider an antibiotic-sensitive (AS) bacterial strain and an antibiotic-resistant (AR) bacterial strain -termed strain 1 and strain 2, respectively-spreading among patients, such that the infection by one bacterial strain provides immunity against the other. Because antibiotics are commonly used in hospitals to prevent a wide range of conditions, Lipsitch et al. [2] assume that patients in the ward are routinely provided antibiotics 1 and 2, regardless of these patients being or not infected by bacteria; more concretely, antibiotic 1 is only effective against the AS bacterial strain, while antibiotic 2 is effective against both strains of bacteria. The acquisition of resistance by bacteria can lead to some fitness cost, amounting to a reduction of the bacterial strain infectiousness due to the corresponding mutation; to represent this fact, Lipsitch et al. [2] consider a common infection rate β = 1.0 days −1 , and set β 1 = β and β 2 = (1 − c)β with c ∈ (0, 1). Spontaneous clearance of sensitive and resistant bacteria occurs at a rate γ, and contributions of antibiotics 1 and 2 to this recovery are represented by rates τ 1 and τ 2 . Patients are assumed to be admitted by and discharged from the hospital ward at a common rate µ.
In our numerical experiments (Tables I-VII), we consider a hospital ward with N = 20 patients, initial numbers (I 1 , I 2 ) = (1, 1) of infectives, and values c ∈ {0.05, 0.1, 0.25} of fitness cost. It should be pointed out that, unlike the paper [2] where the deterministic model is related to frequencies, we shall from now on consider rates β 1 = N −1 β and β 2 = N −1 (1 − c)β, since the random variables in the underlying LD-QBD processes X ([1, Section 3.1]), and X 1 and X 2 ([1, Section 3.2]) amount to numbers of infectives.  Table I. Means and standard deviations of the time T until the end of the epidemic spread, and of the numbers I 1 (T ) and I 2 (T ) of type-1 and type-2 infectives when this occurs in the SI 1 , I 2 epidemic model. (Strain 1: antibiotic-sensitive; strain 2: antibiotic-resistant.) In Tables I-III, the interest is in a preliminary scenario with τ 1 = τ 2 = 0.0 (no usage of antibiotics), γ = 0.0 (no spontaneous recovery) and µ = 0.0 (no arrival or departure of patients during the outbreak), which is readily translated into an SI 1 , I 2 epidemic model. For practical use, it is worth noting that the derivatives of a predetermined descriptor D in Tables I-III (i.e., expected values and standard deviations of T , I 1 (T ) and I 2 (T )) satisfy since rates β 1 and β 2 depend on the value c of fitness cost and the rate β; note that Eqs.
In Table I, we compute the values of the mean length E[T ] of an outbreak and the mean numbers E[I 1 (T )] and E[I 2 (T )] of patients infected by AS and AR bacterial strains, respectively, during the outbreak, together with the corresponding standard deviations. Our interest here is in analyzing the impact that small perturbations in the parameters of (β 1 , β 2 , c, β) have on these summary statistics. Partial derivatives of interest are given in Table II of patients infected by the strain of AR bacteria. • Perturbations in the common infection rate β do not affect the random variables I 1 (T ) and I 2 (T ) at all, since these perturbations lead to equal relative changes in β 1 and β 2 , so that positive and negative effects on these variables are balanced out. This is directly related to the fact that the dynamics of the SI 1 , I 2 epidemic model are governed by the ratio β −1 2 β 1 (in this case, becoming (1 − c) −1 ), and not by the particular magnitudes of β 1 and β 2 . • Stochastic uncertainty -represented by σ(T )decreases with increasing values of β 1 , β 2 and β (roughly speaking, the fastest infections occur, the less volatile the length of the outbreak is), and with decreasing values of c. On the other hand, uncertainty about by σ(I 1 (T ))increases with β 2 (since ∂σ(I 1 (T ))/∂β 2 > 0) and decreases with β 1 and c, for similar reasons; note that analogous comments can be made on the strain of AR bacteria in terms of σ(I 2 (T )).
We stress that comments above refer to the sign of the derivatives in Table II and apply regardless of the particular value of c ∈ {0.05, 0.1, 0.25}. A more detailed comparison between derivatives in absolute terms can only be carried out after computing their dimensionless counterparts, which are related to elasticities (i.e., analyzing (θ −1 D) −1 ∂D/∂θ instead of ∂D/∂θ for any descriptor D and parameter θ). These elasticities are given in Table III, with the following insights: • The mean length of the outbreak is more affected by perturbations in β 1 than in β 2 , and this difference is more significant with increasing values of c. This behavior is directly related to the fact that β 2 < β 1 , since c is strictly positive. In fact, we would expect to obtain a value for the elasticity of E[T ] with respect to β 2 (i.e., /∂β 2 ) equal to its counterpart with respect to β 1 (i.e., Elasticity(E[T ]; β 1 )), in the special case c = 0. Moreover, the expected length of an outbreak is inversely proportional to β -represented by Elasticity(E[T ]; β) = −1-, which is to be expected since β −1 = 1 day, where 1 day amounts to the unit of time, and thus the time unit used for E[T ]. • Some symmetries can be identified; for example, it is seen that for any value c of fitness cost. This is explained again by the fact that dynamics in SI 1 , I 2 epidemic models are governed by the ratio β −1 2 β 1 , so that the mean number of patients suffering infection by the AS bacterial strain can increase either by increasing the value of β 1 or decreasing the value of β 2 ; similar comments apply to the expected number E[I 2 (T )], and standard deviations σ(I 1 (T )) and σ(I 2 (T )). • In general, the rate β represents the most important parameter for the random index T , while β 1 and β 2 are equally important for the random variables I 1 (T ) and I 2 (T ), regardless of the value of c. We now incorporate discharge and recovery of patients into the model of Lipsitch et al. [2] by making use of the SI 1 , I 2 S epidemic model with recovery rates γ 1 = γ + τ 1 + τ 2 + µ and γ 2 = γ + τ 2 + µ, and values τ −1 1 = 5 days and τ −1 2 = 10 days, when discharge of patients -who are replaced by susceptible patients-occurs in average in 7 days (i.e. µ −1 = 7 days), and spontaneous recovery occurs in average in 30 days (i.e., γ −1 = 30 days). These values for τ 1 , τ 2 , µ and γ correspond to realistic selections used by Lipsitch et al. [2, Figure 2], but parameters are known to vary within concrete ranges. For instance, the average duration µ −1 in hospital stay and the average time γ −1 until spontaneously clearance of bacterial carriage vary between 7 and 20 days, and between 30 and 60 days, respectively; see [2, Table 1]. We also select rates λ 1 = N −1 0.1 and λ 2 = N −1 0.1 to represent infections not directly caused by infectious contacts (for example, due to environmental contamination of the hospital ward), but we should point out that these parameters are an addition not considered explicitly in [2]. Note that, once the partial derivatives of a certain descriptor D in Tables IV-VII (i.e., expected values and standard deviations of the key measures T , X max , T (1), X max (1), T (2), X max (2), I 1 (∞) and I 2 (∞)) are obtained with respect to the primary parameters of θ = (β 1 , β 2 , λ 1 , λ 2 , γ 1 , γ 2 ) T , we can then obtain derivatives of D with respect to secondary parameters (i.e., c, β, τ 1 , τ 2 , µ and γ) as with ∂D/∂c and ∂D/∂β verifying Eqs. (1)-(2), and ∂D/∂µ = ∂D/∂γ = ∂D/∂τ 2 . For the sake of brevity, results in Tables IV-VII correspond only to the choice c = 0.25. We note that long outbreaks obtained in our numerical results (lasting for years;  (2)] ∼ 11 patients infected by the AR bacterial strain. Although implementing control measures within the hospital ward would contribute to decrease the values of these summary statistics (more particularly, Lipsitch et al. [2] consider control strategies such as implementing barrier precautions, improving hand washing compliance levels by healthcare workers, or increasing drug dosage when bacteria is detected in the ward), considering such control actions is out of the scope of this paper, and we focus instead on the local sensitivity analysis for the parameters when no intervention is considered.
In Table V we list values of the partial derivatives of summary statistics with respect to the parameters of (β 1 , β 2 , λ 1 , λ 2 , γ 1 , γ 2 ), in the case c = 0.25. Again, the focus in this table is on the sign of these derivatives, while analyzing these quantities in absolute terms requires a previous normalization in terms of elasticities. As the reader may observe, the mean length of the outbreak increases with increasing values of β 1 , β 2 and λ 2 , and with decreasing values of γ 1 and γ 2 , as one might expect. However, it is also seen that ∂E[T ]/∂λ 1 < 0, which suggests that external infections of patients by the strain of AS bacteria act here as a global protection in the hospital ward, reducing the length of the global outbreak. This can be better explained by analyzing scenarios with smaller and larger values of the fitness cost c. For example, for c = 0.1 (results not reported here), we find that ∂E[T ]/∂β 1 and ∂E[T ]/∂λ 1 are strictly negative, so that when the AR bacterial strain is infectious enough, any kind of infection by the strain of AS bacteria acts as a protection measure for the hospital ward in general terms, that is, when analyzing the global outbreak length E[T ]. On the other hand, in the case c = 0.5 -when the fitness cost is large and, as a result, the AR bacterial strain is not so infectious-we find that ∂E[T ]/∂β 1 and ∂E[T ]/∂λ 1 are strictly positive, representing the fact that the protective role of the AS bacterial strain is not worth it here, given the low infectiousness of the AR bacterial strain. The scenario in Table V (c = 0.25) should be considered as an intermediate situation, where external infections by AS bacterial strain help to protect the ward, while infectious contacts among patients by AS bacterial strain do not play the same protective role. These results suggest that, when considering the implementation of control strategies, more focus on avoiding environmental contamination by the strain of AS bacteria, or on avoiding infectious contacts between patients by AS bacterial strain (through healthcare workers), should be made depending on the suspected infectiousness of the strain of AR bacteria also present in the ward. Other insights from Table V are as follows:  (2)], but also to the expected steady-state numbers E[I 1 (∞)] and E[I 2 (∞)] of patients infected by AS and AR bacteria, respectively. • Derivatives of the standard deviation are difficult to interpret in the SI 1 , I 2 S epidemic model. However, we may note here that, for example, the partial derivatives ∂σ(X max )/β 1 , ∂σ(X max )/β 2 , ∂σ(X max )/λ 1 and ∂σ(X max )/λ 2 are strictly negative and, on the contrary, ∂σ(X max )/γ 1 and ∂σ(X max )/γ 2 are strictly positive. This means that the peak of the outbreak behaves in a more deterministic way with increasing values of the infection rates, while it behaves more stochastically when infection and recovery rates are more balanced, that is, with increasing values of recovery rates in Table V.
In order to identify the most important parameters for each descriptor, we list in Table VI the elasticities of summary statistics with respect to the parameters. An examination of Table VI reveals the following observations: • Symmetries identified in the SI 1 , I 2 epidemic model disappear in the multi-type SIS case, since incorporating recovery of patients in the hospital ward into the SI 1 , I 2 S epidemic model results in a more complex description. In particular, the specific magnitudes of β 1 and β 2 have a significant impact on the descriptors, regardless of maintaining the same value for the ratio β −1 2 β 1 . • When analyzing the expected length E[T ] of the global outbreak, the magnitudes of β 2 and γ 2 are the most relevant ones. This fact is closely related to our comment above where the main contribution to the global outbreak length corresponds to long AR bacterial strain outbreaks. The expected length E[T (1)] of AS bacterial strain outbreaks is more affected by β 1 and γ 1 , while β 2 and γ 2 have more impact on the expected length E[T (2)], as  (2)]. • On the other hand, the global peak E[X max ] of infection, as well as the steady-state numbers E[I 1 (∞)] and E[I 2 (∞)] of infected patients by strains of AS and AR bacteria, respectively, seem to be approximately equally affected by infection and recovery rates corresponding to both bacterial strains, which results in comparable absolute magnitudes for elasticities of these descriptors with respect to β 1 , β 2 , γ 1 and γ 2 .
Finally, we list in Table VII the elasticities of descriptors with respect to the parameters of (τ 1 , τ 2 , γ, µ, c, β), with the following insights: • The most relevant parameters correspond to µ, c and β for most of the summary statistics. This means that the discharge of patients (who might be carrying the bacteria), the infectiousness of the bacterial species (represented by β) and the fitness cost c of the antibiotic-bacterial strain are the most important factors affecting the dynamics of these infections. Since discharge of patients carrying the bacteria implies clear ethical concerns not addressed in Ref. [2], and parameters c and β correspond to factors that can not be controlled by policy-makers in the hospital ward, these elasticities are less helpful in terms of analyzing the efficacy of potential control strategies. Policy-making related parameters in this model correspond to the usage of antibiotics 1 and 2.