Sensitivity analysis of the recovery time for a population under the impact of an environmental disturbance

Wildlife populations are often affected by natural or artificial disasters that reduce their vital rates leading to drastic fluctuations in population dynamics. We use a stage‐structured matrix model to study the recovery process of a population given an environmental disturbance. We focus on the time it takes the population to recover to its pre‐event level and develop general formulas to calculate the sensitivity and elasticity of the recovery time to changes in the initial population, vital rates, and event severity. Our results suggest that the recovery time is independent of the initial population size but it is sensitive to the initial structure. Moreover, the recovery time is more sensitive to reductions in vital rates than to the duration of the impact of the event. We explore an application of the model to the sperm whale population in Gulf of Mexico following a disturbance such as the Deepwater Horizon oil spill.


Abstract
Wildlife populations are often affected by natural or artificial disasters that reduce their vital rates leading to drastic fluctuations in population dynamics. We use a stagestructured matrix model to study the recovery process of a population given an environmental disturbance. We focus on the time it takes the population to recover to its pre-event level and develop general formulas to calculate the sensitivity and elasticity of the recovery time to changes in the initial population, vital rates, and event severity. Our results suggest that the recovery time is independent of the initial population size but it is sensitive to the initial structure.
Moreover, the recovery time is more sensitive to reductions in vital rates than to the duration of the impact of the event.
We explore an application of the model to the sperm whale population in Gulf of Mexico following a disturbance such as the Deepwater Horizon oil spill.

Recommendations for Resource Managers
• Understanding a population's recovery process following a disturbance is important for management and conservation decisions.

INTRODUCTION
Over the centuries, the rate at which human activity affects the ecosystem has been accelerating (Vitousek, Mooney, Lubchenco, & Melillo, 1997). The ever increasing human population size and higher demands for resources are reducing biological diversity, causing climate change, and transforming the landscape. In particular, rates of species extinction are now on the order of 100 to 1000 times those before humanity's dominance of Earth (Pimm, Russell, Gittleman, & Brooks, 1995). Direct impacts on wildlife populations come in the form of human-caused disasters (such as oil spills), overexploitation, invasive species introduced by human movement, and habitat fragmentation or degradation (Fahrig, 2017;Hutchings, 2005;Vitousek et al., 1997). In addition, indirect impacts such as shrinking habitat, climate change, and chemical pollution harm species survival in the long term (Frumhoff, 1995;Gibbons et al., 2000;Isaak et al., 2010). The increasing awareness of this situation has inspired more and more conservation efforts to protect wildlife and aid the recovery of threatened species (Brown, 2010). Mathematical modeling is a powerful tool for understanding wildlife population dynamics and setting guidelines for such conservation efforts. The development of mathematical models that can be parameterized by empirically derived data is particularly important. Numerous population models have been developed to explain and predict population dynamics and explore efficient management strategies for at risk species (Adams, Banks, Banks, & Stark, 2005;Baveco, & De Roos, 1996;Cushing, & Zhou, 1994). Among these, matrix population models are very popular with both mathematicians and biologists. This is not only because of their simplicity and ability to directly link vital rates with population dynamics (Fieberg & Ellner, 2001), but also because there exists well-established theories of perturbation analysis and stability analysis using matrix population models (Caswell, 2001(Caswell, , 2012Shyu & Caswell, 2014). For instance, Wisdom and Mills (1997) used matrix population models to study the declining greater prairie-chickens population. Applying sensitivity analysis, they concluded that the combination of nest success and brood survival is the most effective way to increase the growth rate of the population. In another study, Bridges, and Carroll (2000) used a matrix population model to project the ecological effects of chronic sediment toxicity on the estuarine amphipod Leptocheirus plumulosus. Both studies assumed that the environment was constant and focused on examining the asymptotic growth rate and its sensitivity to model parameters, which provides useful information such as the population trend and reveals the most important life stages of Natural Resource Modeling 3 of 22 the population. However, such analysis cannot be applied to environment dependent matrix population models in which the projection matrices change over time. In this paper, we introduce a new population endpoint, the recovery time, which can be useful in such situations. For a slow growing population, a small environmental disturbance may result in a declining population (Chiquet, Ma, Ackleh, Pal, & Sidorovskaia, 2013). In this paper, we examine the recovery ability, the time it takes for such population to recover, and how different factors contribute to the recovery process. Understanding a population's recovery process following a disturbance can aid in management decisions, such as establishing regulations on harvesting or habitat quality, and may also help in predicting what may happen to a population should another disturbance occur. There are various measures for population recovery including the species reappearance, the return to predisturbance total density or biomass, and the average size of individuals in a species (Detenbeck, DeVore, Niemi, & Lima, 1992). Here, we define the recovery time as the time it takes for an affected population to reach its predisturbance level. Naturally, there are many factors affecting the recovery time. The vital rates during the recovery process are particularly critical to the length of recovery time since the impact of many other factors is realized through them. By estimating the recovery time for different values of vital rates, one can better understand the recovery process and subsequently choose reasonable intervention strategies. We derive formulas for the sensitivity of the recovery time to changes in the initial population structure, vital rates, and environmental factors. This sensitivity analysis pinpoints which factors contribute the most to the recovery process. The techniques used in the sensitivity derivations are general enough to be applied to various models. They can also be used to obtain the mean recovery time and upper and lower estimates of the recovery time for a population experiencing demographic stochasticity, provided that the recovery probability is equal to or near one. Thus, these formulas provide an alternative to stochastic simulations that has the advantage of being computationally efficient and does not require solving the stochastic model thousands of times to generate a mean recovery time.
This paper is organized as follows. In Section 2, we present a model for a structured population whose vital rates are proportionally reduced for a specified amount of time due to an environment disturbance, after which they return to their pre-event levels. In Section 3, we establish formulas to calculate the sensitivity (elasticity) of the recovery time to the initial population structure, the vital rates, and the environmental parameters. We also derive the formulas to calculate the mean and variance of the population with demographic stochasticity. In Section 4, we apply our results to consider the potential effect of the Deepwater Horizon (DWH) oil spill on the sperm whale population in Gulf of Mexico. In addition, we consider the effect of demographic stochasticity on the recovery time and calculate the confidence interval of the recovery time and its sensitivity or elasticity with respect to the environmental parameters under different scenarios. Finally, in Section 5, we summarize our results and discuss their strength and limitations.

MODEL FORMULATION
We consider a population described by a discrete-time, stage-structured matrix model. We divide the female population at time into stages described by the vector ( ) ∶= [ 1 ( ), 2 ( ), … , ( )] ⊺ , where ⊺ denotes the transpose of a vector. Given the population at time , the population at time + 1 is determined by the projection matrix [ ( ( )), ( )]. Here, ( ) describes the environment at time and ( ( )) represents the environment-dependent vital rates. This leads to the nonautonomous matrix model where = 0, 1, 2, …. Iterating Equation (1), we have ( + 1) = ∏ =0 [ ( ( )), ( )] (0). We set the time when an environmental disturbance happens to be = 0. Following a disturbance, we assume that the population experiences reductions in vital rates. These vital rates are allowed to recover over time as the impact of the disturbance vanishes. We assume that the vital rates prior to the disturbance yield a growing population, since the population will go extinct otherwise even without an external disturbance. We also assume that these reductions are sufficiently large to cause a declining population. We focus on the transient population dynamics and consider the simple case where no additional catastrophic event happens during the recovery process.
A1. Assume the dominant eigenvalue of [ ( ( )), ( )] is greater than one in the absence of a disturbance and is less than one at (0).
In general, the recovery of the vital rates approximately follows a sigmoid curve, but as shown in Ackleh et al. (2017), a constant reduction followed by the complete recovery of the vital rates provides lower and upper bound estimates for such sigmoid functions. Thus, to model the process of environment recovery after an event, we consider the following step recovery function: where 0 represents the proportion of reduction in the vital rates caused by the event. For instance, 0 could represent reductions caused by abnormally high toxicant levels due to an oil spill or low food sources due to flooding. In general, is the time when the negative effect caused by the environmental event disappears.
To examine the recovery process, we focus on the recovery time which we define to be the first time after at which the population has recovered to its pre-event size. From Equations (1) and (2), the matrix model can be written as where 0 represents the pre-event projection matrix (where ( ) = 0) and 0 represents the projection matrix for the impacted population (where ( ) = 0 ). When the environment is described by Equation (2), assumption A1 is equivalent to the following assumption. A2. Assume [ 0 ] > 1 and [ 0 ] < 1, where [⋅] denotes the spectral radius. Using Equation (3), for > , ( + 1) can be written as Let ( ) ∶= ⊺ ( ) denote the total population size at time , where is an × 1 vector of ones. Then, the recovery time, denoted by , is the first integer + 1 > that satisfies ( + 1) ≥ (0) or, from Equation (4), We define the × 1 vector ( ) to be the population distribution of individuals from all stages at time so that ( ) = ( ) ( ). For notational simplicity, we use and to denote the total population and population distribution at time , respectively. Using this notation, Equation (5) is equivalent to where 0 is the population distribution at time = 0. The assumption that [ 0 ] > 1 guarantees that a solution to Equation (6) exists. We are interested in examining how the recovery time is affected by changes in properties of the environment or the population. Two types of perturbation analysis can be used to consider how an outcome changes if a parameter is slightly changed (Caswell & Gassen, 2015). One type is sensitivity analysis, which describes how additive changes in affect when other parameters are held constant and is given by . The other type is elasticity analysis, which describes how proportional changes in affect while keeping other parameters unchanged and is given by . However, both of these analyses require taking the derivative of the discrete variable . To address this issue, we define an auxiliary function of a continuous variable by The function is a real-valued function of ∈ ℝ + whenever 0 has no nonpositive real eigenvalues. In addition, this function is equal to the left-hand side of (6) for all ∈ ℕ. Therefore, in order to understand the sensitivity of the recovery time to various parameters, we can study the auxiliary recovery equation 1 = ( ) instead, that is Proof of the existence of a real solution to Equation (8) whenever 0 has no nonpositive real eigenvalues is provided in the Appendix. The recovery time , as defined by Equation (6), can be obtained by solving Equation (8) and rounding the solution up to the next integer value. Meanwhile, the sensitivity of with respect to a perturbation in a variable , , can be obtained by implicitly differentiating Equation (8), solving for , and evaluating this quantity at a solution pair ( , ) of Equation (8). Of course, in order to interpret these values in the context of the discrete equation (6), these sensitivities should be rounded up to the next integer value. In this manner, the auxiliary equation (8) can be used to investigate how perturbations in parameters affect the solution to Equation (5).
The assumption of no nonpositive real eigenvalues commonly arises when trying to convert a system from discrete to continuous time (see, for instance, Singer & Spilerman, 1976). The process of solving for the recovery time in Equation (8) (as seen in the Appendix) involves taking the logarithm of 0 . This is complex if 0 has any real negative eigenvalues and results in a complex-valued recovery time. Therefore, to ensure that the output is biologically reasonable, we only consider the case where matrix 0 has no nonpositive real eigenvalues. In the following section, we analyze the sensitivity of the recovery time with respect to changes in the initial population, vital rates, and the environment.

Sensitivity of the recovery time to the initial population
In this section, we first investigate the sensitivity of the recovery time to the initial population. In particular, we examine the sensitivity of with respect to the initial population size 0 and the initial population distribution 0 . We show that the recovery time is independent of the initial population size.
Proposition 3.1. Assume A2. The recovery time is independent of initial population size (0).
Since the recovery time is determined by the first ∈ ℕ that satisfies (6), Proposition 3.1 follows from the fact that Equation (6) is independent of 0 . We note that Proposition 3.1 is not specific to the environment as described by Equation (2). When the projection matrices are diagonalizable, the right-hand side of Equation (6) can be decomposed in terms of the eigenvalues and eigenvectors of 0 and 0 , which reduces computation time. Details on this decomposition are provided in the Appendix.
Next, we study how changes in the initial distribution 0 affect the recovery time of the population. To do this, we have to take into account the relationship of the components in 0 . Since 0 is a probability distribution (and thus its components must add to one), if we perturb a single component of 0 , we must compensate by perturbing other components of 0 so that 0 remains a probability distribution. Therefore, we define the derivative of the recovery time with respect to a single component 0 as where 0 0 reflects the compensation strategy. As pointed out in Caswell (2013), there are different ways to define the relation between the components of the distribution vector. In Section 4.1, we consider two particular strategies to define how the components relate to each other. We give the general formula for the sensitivity for an arbitrary relation in the following proposition.
Proposition 3.2. Assume A2 and that 0 has no nonpositive real eigenvalues. The sensitivity of with respect to 0 is given by where is a scalar given by Here, the ⊗ operator denotes the Kronecker product and the vec operator stacks the columns of a matrix to make a column vector.
Proof. We begin by taking the derivative of Equation (8)  Given two functions ∶ ℝ → ℝ × and ∶ ℝ → ℝ × , the chain rule for matrix derivatives (Magnus & Neudecker, 1988) For the derivative of the right-hand side of Equation (8), applying Equation (11) repeatedly, we have Theorem 1.31 in Higham (2008) states that matrices whose eigenvalues do not include nonpositive real values have a unique principal logarithm. Thus, there exists a unique matrix such that 0 = , where is the principal logarithm of 0 whose eigenvalues lie in the strip This theorem allows us to write the term Substituting Equation (13) into Equation (12) yields ( Finally, combining 1 0 = 0 and Equation (14), we have , we solve this equation for and evaluate the result at = to obtain formula (10). ■

Sensitivity of the recovery time to vital rates
From model equations (1) and (2), we see that the recovery time depends on the initial population condition, the environmental parameters, and the vital rates. For endangered species with declining populations, there are various management strategies, such as limiting human hunting or fishing, and predator control efforts that can increase their vital rates (Silvy, 2012). Here, we focus on how improvements in vital rates as a result of management may affect the recovery time. We assume that management efforts are only applied while the vital rates of a species are being negatively impacted by a disturbance and are removed after the vital rates recover to predisturbance values. Let̂∶= [ 1 (0, ), 2 (0, ), … , (0, )] denote a vector of the vital rates impacted by the disturbance, where (0, ) represents a constant vital rate from year zero to year . Then,̂gives the sensitivity of the recovery time to changes in the affected vital rates from year zero to year . This sensitivity analysis reveals which life stages are particularly important for the recovery of the population.
Proposition 3.3. Assume A2 and that 0 has no nonpositive real eigenvalues. The sensitivity of with respect to changes in vital rates up to year is given bŷ Proof. We start by taking the derivative of Equation (8) with respect tô. On the left-hand side, the derivative is zero. On the right-hand side, applying Equation (11), we have For the first term on the right-hand side of Equation (16), we have For the second term on the right-hand side of Equation (16), we apply Equation (11), the chain rule, and the property for a square matrix from Magnus and Neudecker (1988), where is an integer. We obtain . Natural Resource Modeling 9 of 22 Substituting these two terms back into Equation (16) yields Solving for̂and evaluating the result at = , we obtain Equation (15). ■

Sensitivity of recovery time to environmental parameters
Numerous studies of environmental disturbances have shown that the vital rates of wildlife populations in the affected area experience declines following the event (Lane et al., 2015;Mensah, Palmer, & Muller, 2014;Relyea, 2009). However, for some species such as sperm whales and beaked whales in Gulf of Mexico, the habitat and the behavior of these large marine mammals make it difficult to quantify the changes in the vital rates caused by such events or the duration of impact. In other words, it is complicated, sometimes impossible, to obtain accurate estimates of 0 or . Thus, it is useful to know how the recovery time changes if the parameters describing the impact of a disturbance change.
In addition, even if 0 and can be estimated accurately, the study of the sensitivity of the recovery time with respect to the environmental parameters reveals to us how much the recovery time can be reduced if we are able to reduce the impact of the event. Thus, we study how the change in the impact magnitude 0 and duration of impact affect the recovery time .
Proposition 3.4. Assume A2 and that 0 has no nonpositive real eigenvalues. The sensitivity of to changes in 0 is given by Proof. To obtain the sensitivity of the recovery time to 0 , we take the derivative of Equation (8) with respect to 0 . The derivative of the left-hand side of Equation (8) is zero, while the derivative of the right-hand side, obtained by applying Equation (11) repeatedly, is ( and evaluating the result at = , we obtain Equation (17). ■ We also study the sensitivity of with respect to changes in . Since is a discrete parameter, we use the following difference operator: to estimate this sensitivity. Equation (18) gives the average rate of change of given a one unit change in . In Section 4.1, we use formulas (17) and (18) to examine a sperm whale population subject to an environmental disturbance.

The mean and variance of a population with demographic stochasticity
So far, we have only considered the recovery time for a deterministic population model. However, demographic stochasticity, which takes into consideration the variation among individuals in the random outcomes of the vital rates, can also affect the recovery times. When demographic stochasticity is incorporated into a population model, each individual experiences an independent realization of the vital rates, making the population size at any time a random variable. As a consequence, the recovery time is also a random variable that may take different values for each simulation run. To obtain the mean and variance of the population at certain times without performing thousands of simulations, the analysis of the stochastic process, which is based on a multitype branching process, is needed. See Caswell (2001, Chapter 15) for a detailed discussion of the relation between multitype branching processes and matrix population models. We impose the same assumptions on the stochastic demographic events as in Ackleh et al. (2017). We begin by decomposing the projection matrix ( ) into a transition matrix ( ) = ( ( )) and a fertility matrix = ( ( )) such that ( ) = ( ) + ( ). The fertility matrix ( ) is zero elsewhere besides the nonzero birth entries from ( ). The matrix ( ) does not include death as a state, so we add an additional row to ( ) to creatẽ( ) whose columns add up to one.
Let ( ) ( ) ∶= [ ( ) 1 ( ), … , ( ) ( )] ⊺ be the vector random variable that gives the number of offspring of each stage from a parent in stage at time , where represents the number of stages in the population model and ( ) ( ) denotes the random number of offspring of stage produced by a parent of stage at time . Here, we use offspring to refer to both individuals produced through transition and through birth. Then, we can write where the random variables ( ) ( ) and ( ) ( ) give the number of offspring of stage produced by a parent of stage via transition and birth at time , respectively. We assume that the random variable ( ) ( ) follows a multinomial distribution determined by the th column of̃and ( ) follows a binomial distribution determined by the first row and th column of . Given the independence of birth and transition events, the expected number of offspring of stage produced by parents of stage can be obtained by ( ( ) ( )) = ( ( ) ( )) + ( ( ) ( )) = ( ) + ( ) = ( ), where ( ) is the ( , )th component of ( ). We further assume that the births of different offspring stages are independent. Then, we have where Var( ( ) ( )) denotes the variance of the random variable ( ) ( ) and Cov( ( ) ( ), ( ) ( )) denotes the covariance between random variables ( ) ( ) and ( ) ( ). We know that where ( ) represents the vector including offspring of all stages from the th parent of stage . Thus, ( ) are identical random variables. Applying the rule of total expectation and Equation (19), we can write the expected population at time + 1, ( ( + 1)), as We proceed to calculate the variance of the population ( ). We use ( ( )) = ( ( )) to denote the covariance matrix at time for ( ), where ( ) represents the covariance between ( ) and ( ). Using the rule of total expectation and Equation (21), we have ( ( + 1)) = (Cov( ( + 1)| ( ))) + Cov ( ( ( + 1)| ( ))) where the covariance matrix of ( ) can be obtained using Equation (20). Therefore, we have where ( ) is obtained from Equation (22).

APPLICATION TO A SPERM WHALE MODEL
In this section, we apply the results from Section 3 to study the sperm whale population in the Gulf of Mexico following a disturbance such as the DWH oil spill. Sperm whales are known to have been 0.4920 present in areas of the Gulf impacted by the oil spill (Ackleh et al., 2012) and studies have shown that the DWH oil spill has affected the vital rates of large marine mammals in Northern Gulf of Mexico both through the toxins released during the oil spill and the chemical dispersants thrown in (Geraci, 1990;Lane et al., 2015;Wise, Wise, Wise, Thompson, & Wise, 2014a;Wise et al., 2014b). In addition, it has been shown that the chemicals released from the oil spill continue to affect the marine mammals after contact (Geraci, 1990;Wise et al., 2014a, b). We use the results from Sections 3.1-3.3 to study the sensitivity of the recovery time of the sperm whale population following a disturbance. In addition, we use results from Section 3.4 to investigate how the recovery time is affected by demographic stochasticity and obtain an interval estimate of the recovery time under different scenarios.

Sensitivity analysis of recovery time
Following Chiquet et al. (2013), we divide the sperm whale population into five stages: calves, juveniles, mature females, mothers, and postbreeding females. For notational simplicity, we use ( ) ∶= [ ( ( ))] to denote the projection matrix at time where ( ) = ( )(1 − ( )) represents the probability of surviving and staying in stage at time , = ( ) ( ) represents the probability of surviving and moving to stage + 1 for = 1, … , 4, and 5 = 5 ( ) 5 ( ) represents the probability that a postbreeding female survives and returns to the mature female stage. Here, ( ) represents the survival probability of whales in stage at time and ( ) represents the probability of whales in stage moving out of stage at time . The mature female whales have fecundity = 0.125. The values of the survival rates and transition rates at their pre-event levels are listed in Table 1.
Here, we consider the case where a spill results in a reduction in survival rates so that ( ( )) = [ 1 ( ( )), … , 5 ( ( ))], where ( ) is defined in Equation (2). Given that juvenile stages are known to be more sensitive to toxicants in many species (Birge, Black, & Westerman, 1979), we allow for the possibility that a disturbance has a greater impact on the survival of the juvenile stages. We define the reduction in the survival of stage according to Since Ackleh et al. (2017) observed that the recovery probability for a sperm whale population is not heavily impacted by these weights, for simulation purposes, we consider only one scenario with = 2 for = 1, 2 and = 1 for = 3, 4, 5.
Using the parameters above, we first calculate the recovery time for different population sizes when the environment is given by = 10 and 0 = 0.05. As expected from Proposition 3.1, the recovery times for all the simulations are the same (not shown). Next, we consider the sensitivity of the recovery time of the sperm whale population to perturbations in the initial population distribution. Assume that the female to male ratio of sperm whales in Gulf of Mexico is one to one. Choosing an initial female population size of 0 = 381 (Waring, Josephson, Maze-Foley, & Rosel, 2012), we start by considering how the total population size at time , ( ), changes with respect to different initial distributions. Figure 1 shows the total population size for = 1, … , 100 with six different initial distributions 0 . Here, we define 0 = for = 1, … , 5, where is a vector with 1 on the th entry and zeros elsewhere, and we define 6 0 to be the stable stage distribution 6 0 = [0.0850 0.2077 0.3617 0.1783 0.1672] ⊺ . From Figure 1, we see that the total population with the initial distribution 3 0 is always larger than the others. In particular, there is 170% more population with 3 0 than 1 0 at = 100. This implies that, given a fixed pre-event population level, the population can recover in a shorter time if a larger portion of the population consists of mature female adults. We note that, since the sperm whale model is linear, it is possible to increase or decrease the number of stages in the model while preserving the dominant eigenvalue and stable stage distribution (Salguero-Gomez & Plotkin, 2010). Therefore, the long-term dynamics of the population are not impacted by the dimension of the model. However, since the recovery time depends on the transient dynamics of the model, as can be seen in Figure 1, changing the dimensions of the model may have an impact on the recovery time.
To better understand how the initial population distribution affects the recovery time, we apply Proposition 3.2 to the sperm whale population subject to reductions defined by (25). We consider the following two strategies to specify the relation between the components of 0 . In Strategy 1, we distribute the change in 0 uniformly among the other components, in which case 0 0 = − 1 4 for ≠ (Caswell, 2001). In Strategy 2, we create an × matrix where the ( , ) component represents the sensitivity of with respect to 0 given that the change in 0 is compensated by the change in 0 . In  Figure 2, we see that, using Strategy 1, increasing the proportion in either of the two juvenile stages leads to a longer recovery time, while increasing the proportion in any of the three adult stages shortens the recovery time. This agrees with what we observe in Figure 1 in the sense that when there are more individuals in the adults stages, it takes less time to recover to a fixed population level.
Overall, we see that the mature female adult population reduces the recovery time the most. The right side of Figure 2 gives the elasticities of the recovery time when we apply Strategy 2 to the stable stage distribution. The entry in the th row and th column describes how the recovery time changes when an increase in stage is compensated by a reduction in stage . We conclude that increasing the stage 3 proportion and decreasing the stage 1 proportion has the greatest effect on the recovery, reducing the recovery proportion by 36%.
To understand how the vital rates affect the recovery time under different magnitudes of reduction and impact time, we use Proposition 3.3 to compute the elasticity of the recovery time to the survival rates. Figure 3 gives the elasticities of the recovery time with respect to the survival rates for two different values of . We observe that these elasticities are negative. This is biologically reasonable since higher survival rates mean larger growth rates, which result in the population reaching the preevent level faster. In addition, our sensitivity calculations show that as 0 increases, an increase in any of the three mature stages yields a larger reduction in recovery time for larger 0 . Meanwhile, an increase in either of the two immature stages leads to a smaller reduction in recovery time for larger 0 . For the succinctness of the presentation, these figures are not shown. However, from Figure 3, we see that, as 0 increases, an increase in the survival rate of any of the five stages leads to a smaller proportional reduction in recovery time. This is because even though the absolute change in reduction time is bigger, the proportion of change is actually smaller.

Examining the effect of demographic stochasticity on the recovery process
Utilizing the mean and variance formulas given in Equations (21) and (22), we calculate lower and upper estimates for the recovery time, and ℎ ℎ , from the 95% confidence interval of the total population size. In Dennis and Patil (1984), it is shown that, for many stochastic population models, the gamma distribution arises as a distribution of population size when the equilibrium population Natural structure is reached. By assuming that each stage of the population follows a gamma distribution with parameters obtained using the mean and variance from Equations (21) and (22), it is possible to analyze the sensitivity of recovery time using the stochastic model. We perform the following procedure to obtain and ℎ ℎ for different initial population sizes.
Simulation Procedure: (i) Initialize parameters in model (3)  (iii) Use Equations (21) and (23) to calculate the mean and variance of each stage at time = . Assuming the population follows a gamma distribution, obtain two vectors, 0 and ℎ ℎ 0 , whose entries contain the 95% confidence intervals for each stage at time = .
(iv) Using the map ( + 1) = 0 ( ), find the first time the total population size reaches 0 using initial conditions 0 and ℎ ℎ 0 . Add to the solution with initial condition 0 to obtain ℎ ℎ . Add to the solution with initial condition ℎ ℎ 0 to obtain .
From Figure 4, we see that, given a fixed initial population distribution, the recovery time calculated from the 95% confidence interval of the population size becomes less sensitive to the initial population size as the initial population size increases.
We use the Simulation Procedure to obtain the recovery time calculated from the 95% confidence interval of the population size for different reduction proportions 0 by changing step (v) to "Repeat step (i) to step (iv) for 0 = 0.01 ∶ 0.02 ∶ 0.30," and using initial female population 0 = 381. As we can see from the left side of Figure  increases. This agrees with our expectation since, as the mean recovery time increases, the variability during the recovery period grows. Note that the mean recovery time, , and ℎ ℎ all increase almost linearly with respect to 0 for 0 < 0.2. The right side of Figure 5 gives the elasticity of the recovery time with respect to 0 . From the right side of Figure 5, we see that the elasticity of the recovery time increases when 0 increases. This tells us that the recovery time grows superlinearly when 0 increases.
We also investigate how affects the recovery time and the interval [ , ℎ ℎ ] by changing step (v) to "Repeat step (i) to step (iv) for = 2 ∶ 2 ∶ 40," and using initial female population 0 = 381. From the left side of Figure 6, we observe that as increases, the recovery time increases and both and ℎ ℎ become larger. The right side of Figure 6 gives the sensitivity of the recovery time with respect to changes in calculated using Equation (18). We see that the recovery time increases linearly with a factor of 6 when increases. In other words, it takes six more years to recover to the pre-event population level if the survival rates are reduced for an additional year. Compared to the effect of changes in the impact magnitude 0 on the recovery time, the effect of changes in on the recovery time is very small.
We note that, if the environmental parameters are such that there is a high probability a stochastic population will go extinct, then simulation procedure may not be appropriate for examining population recovery. The reason for this is that the implicit formula (8) assumes that populations always recover.
In such a scenario, it may be more appropriate to use stochastic simulations to calculate the recovery time where, for populations that never recover, the recovery time is set at an arbitrarily large value. This would result in a mean recovery time and interval [ , ℎ ℎ ] that grow much faster than those calculated using Equation (8). For comparison purposes, we also calculated the 95% confidence interval for the recovery time using a simulation procedure similar to the one described in Ackleh et al. (2017). For the range of parameter values considered in Figures 4-6, and ℎ ℎ closely match this confidence interval, with values within 1% of each other (not shown).

DISCUSSION AND CONCLUSION
In this paper, we present a stage-structured nonautonomous matrix population model and use it to study the recovery process of a vulnerable population going through some environmental disturbance, such as overexploitation, an oil spill, or chemical pollution. We assume that the disturbance causes a reduction in vital rates that results in a declining population size. After the effect of the disturbance has dissipated, we assume that the vital rates return to predisturbance values. Using a step function to describe this change in vital rates, we establish an implicit equation for the recovery time given in Equation (5) and examine the sensitivities of the recovery time using the auxiliary equation (8). Using methods from matrix calculus, we derive formulas for the sensitivity of the recovery time with respect to changes in the vital rates, initial population condition, and environmental factors. In addition, utilizing Equations (21) and (23), we estimate upper and lower bounds of the mean recovery time for a population experiencing demographic stochasticity. Our analysis relies on the fact that the projection matrix 0 has no nonpositive real eigenvalues and it remains an open question whether similar calculations will hold when this assumption is removed.
As an example, we demonstrate how the sensitivity formulas can be used to study the sperm whale population in Gulf of Mexico following a disturbance such as the 2010 DWH oil spill. As shown in Figure 4, a small population is associated with a large variance in the mean recovery time. The reason for this is that demographic stochasticity has a greater influence on small population, increasing the possibility of significant differences for different realizations. Since the sperm whale population in the northern Gulf of Mexico is relatively small, which means that demographic stochasticity is likely to influence the recovery time, we calculate the recovery time from the 95% confidence interval of the total population size for a range of environmental conditions, defining an interval [ , ℎ ℎ ]. Figures 5 and 6 show that there is a large variation in the recovery time and that the interval [ , ℎ ℎ ] becomes wider for more severe reductions in vital rates. In this application, the limited information regarding the vital rates are obtained from Chiquet et al. (2013). The stock assessment of the sperm whales in Northern GoM from the literature does not include information about the structure of their population. Even less is known about the quantitative impact of the oil spill on the sperm whale population. Our calculation of the elasticity of recovery time to the vital rates shows that, over a wide range of values, uncertain or erroneous estimates of just a few rates would have a serious impact on the recovery time. We find that the mature stages, especially the reproducing female stage, contribute the most to shortening the recovery time. Comparisons of Figures 2 and 3 reveal that the initial population condition does not affect the recovery time nearly as much as the vital rates.
The recovery time analysis presented in this paper assumes that the disturbance causes a constant reduction in the vital rates for years, after which the effect disappears. Even though this is an oversimplification, by adjusting , the calculations presented in this paper can provide lower estimates or upper estimates for a more realistic scenario. Thus, the recovery time can be used by management organizations as an estimation for conservation efforts. For example, after years of overexploitation of marine fishes, limited fishing can be enforced until the upper estimate of the recovery time to allow the stock to recover. In addition, since the lengths of marine fishes provide good estimates of their ages (Goodyear, 1995), the sensitivity analysis of the recovery time can provide suggestions for the restrictions on sizes of marine fishes for commercial fishing.
In other cases where a disturbance is caused by an accident such as an oil spill, the assumption that additional accidents do not occur before the population is recovered may not hold. However, the recovery time defined in this paper may still be used as the best scenario estimation since any further disturbance may lead to reduced vital rates based on the existing level, which could slow down or even reverse the recovery process. Furthermore, from results in Section 4, we expect the sensitivity analysis to reveal the same significant stages even if more accidents occur. If more disturbances occur before the population is able to recover, we may expect a higher extinction probability. In such a case, the recovery time calculated from Equation (8) is no longer a reliable approximation to the recovery time obtained from the stochastic simulations. It is in our interest to develop creative methods and techniques in future studies to address this problem. Proof. Consider the auxiliary function defined by Equation (7). An arbitrary matrix power is defined via the Cauchy integral (Higham, 2008)  Proof. Since the matrices 0 and 0 are diagonalizable, we write them as where and 0 are diagonal matrices containing the eigenvalues of 0 and , 0 of 0 on their diagonal, respectively. The columns of matrices and 0 contain the corresponding right eigenvectors, = ( 1 2 ⋯ ), 0 = ( 1, 0 2, 0 ⋯ , 0 ), and the rows of −1 and −1 0 contain the corresponding left eigenvectors, Using Equation (A3), Equation (8) can be rewritten as , 0 , 0 , 0 0 , which gives the desired results. ■