Noise source importance in linear stochastic models of biological systems that grow, shrink, wander, or persist

While noise is an important factor in biology, biological processes often involve multiple noise sources, whose relative importance can be unclear. Here we develop tools that quantify the importance of noise sources in a network based on their contributions to variability in a quantity of interest. We generalize the edge importance measures proposed by Schmidt and Thomas [1] for first-order reaction networks whose steady-state variance is a linear combination of variance produced by each directed edge. We show that the same additive property extends to a general family of stochastic processes subject to a set of linearity assumptions, whether in discrete or continuous state or time. Our analysis applies to both expanding and contracting populations, as well as populations obeying a martingale (“wandering”) at long times. We show that the original Schmidt-Thomas edge importance measure is a special case of our more general measure, and is recovered when the model satisfies a conservation constraint (“persists”). In the growing and wandering cases we show that the choice of observables (measurements) used to monitor the process does not influence which noise sources are important at long times. In contrast, in the shrinking or persisting case, which noise sources are important depends on what is measured. We also generalize our measures to admit models with affine moment update equations, which admit additional limiting scenarios, and arise naturally after linearization. We illustrate our results using examples from cell biology and ecology: (i) a model for the dynamics of the inositol trisphospate receptor, (ii) a model for an endangered population of white-tailed eagles, and (iii) a model for wood frog dispersal. Author summary Biological processes are frequently subject to an ensemble of independent noise sources. Noise sources produce fluctuations that propagate through the system, driving fluctuations in quantities of interest such as population size or ion channel configuration. We introduce a measure that quantifies how much variability each noise source contributes to any given quantity of interest. Using these methods, we identify which binding events contribute significantly to fluctuations in the state of a molecular signalling channel, which life history events contribute the most variability to an eagle population before and after a successful conservation effort rescued the population from the brink of extinction, and which dispersal events, at what times, matter most to variability in the recolonization of a series of ponds by wood frogs after a drought.

Noise source importance evaluates the relative contributions of a set of noise sources to 2 variability in a quantity of interest. By evaluating noise source importance we can 3 quantify when, where, and to which quantities, individual noise sources matter. 4 Noise source importance was originally introduced to identify unimportant noise 5 sources that could be ignored as a form of model reduction [1,2]. In some cases, such as 6 in the Hodgkin-Huxley model for neuron firing, this approach can lead to efficient, yet 7 accurate, approximate simulation [3][4][5]. To illustrate this application we will ask, which 8 binding events can be ignored when simulating the inositol triphosphate channel 9 responsible for calcium induced calcium release? 10 More broadly, when fluctuations in a quantity are important, it is helpful to know 11 where those fluctuations come from. For example, variance in the reproductive success 12 of individual organisms can be decomposed into variability contributed by individual 13 traits, and variability arising from random events during the life of the individual. The 14 relative importance of these two sources determine whether an individual's success is 15 determined primarily by their intrinsic quality (pluck), or luck during their lifetime. 16 Snyder et al. applied this decomposition to show the lifetime reproductive success of a 17 shrub and perrenial grass (Artemisia tripartita and Pseudoroegneria spicata) are more 18 determined by luck than the quality of their location [6]. 19 The same kinds of questions can be asked for other noise sources and quantities of 20 interest. In most conservation efforts, population size is the most important quantity. 21 Variability in population size drives stochastic extinction for small populations, so 22 controlling variability is important for their conservation. As an example, we will ask, 23 did variation in clutch size, juvenile survival, or adult survival play a more important 24 role in variability in a population of white-tailed eagles (Haliaeetus albicilla) while they 25 were endangered? Have the important sources of variability changed now that the eagle 26 populations have recovered? Fluctuations also play an important role in colonization 27 and invasion processes triggered by rare dispersal events. We will ask, does variability 28 due to dispersal play an important role in wood frog (Rana sylvatica) populations when 29 at equilibrium, and during a recolonization process following a drought? approximations may not be meaningful if the process is expected to leave the 98 region where the linear approximation is valid. Then it is useful to compute noise 99 source importance at intermediate times. Moreover, in some cases we are 100 interested in the variability along a transient trajectory. For example, when 101 considering an invasion process, we are typically interested in dynamics during the 102 transient approach to equilibrium, not at equilibrium. 103 The paper is organized as follows. In §Scope we outline the sufficient conditions 104 required for a model to suite the analysis. These conditions define the scope of models 105 considered. We then present some important model classes that satisfy the sufficient 106 conditions. 107 We divide our subsequent results into theory and application. In §Theoretical 108 Results we derive the noise source importance measures for discrete-and continuous-109 time models, at finite and long times, in each possible limiting scenario, and show that 110 Schmidt and Thomas' measure [1] is a special case of our generalized measure. We 111 demonstrate that the measures all share a generic derivation in terms of the production 112 and propagation of variance, and the distinctions between the measures arise from 113 whether most present variability was produced recently, or in the far past. We then 114 compute the long time limit of the important measures in each limiting scenario, and 115 find that, in contrast to the existing literature, when the process grows or wanders, long 116 time noise source importance is independent of the choice of observable (c.f. Lemma ). 117 We conclude our theoretical investigation by presenting the equivalent measures in the 118 affine case. 119 In §Applications we apply our theoretical results to illustrative examples 120 corresponding to sample model types introduced in §Scope. The applications illustrate 121 how the noise importance measures can be computed for widely used classes of 122 stochastic models. We consider applications to cellular signalling to demonstrate 123 consistency with existing methods in a model system where noise source importance has 124 not been computed before. A reader interested in efficient simulation and model 125 reduction can skip to §Inositol Trisphosphate and Calcium Signalling on page 15. We 126 also consider applications to white-tailed eagles to demonstrate that our methods can 127 be applied to discrete-time matrix models without a steady state population, and to a 128 nonlinear model of wood frog dispersal to illustrate how our methods can be applied to 129 linear approximations both about equilibrium and along transients associated with Let X(t) denote a continuous-time stochastic process and X t a discrete-time stochastic 136 process. We require the following properties: 137 P1 Markov: The stochastic process is Markovian.

138
P2 Linearity of Conditional Expectation: If X(t) is a continuous-time model, 139 then d dt E[X(t)], given X(t) = x, is either an affine or linear function of x. If X t is 140 a discrete-time process, then E[X t+∆t ], given X t = x, is an affine or linear can be expressed as a sum over noise sources, the derivative is well-defined 148 (symmetric and positive-definite) for any partial sum over the noise sources, and 149 the derivative is zero if no noise sources are included. If X t is a discrete-time 150 process then Cov[X t+∆t ] can be expressed as a sum over noise sources, the 151 covariance is well defined (symmetric and positive-definite) for any partial sum 152 over noise sources, and the covariance is zero if no noise sources are included. In 153 either case the contribution from each noise source is an affine function of x.

154
The first three conditions ensure that the state covariance at time t takes a canonical 155 form. The last condition ensures that the state covariance can be separated into 156 contributions from each noise source. The results provided in the main text assume that 157 both the conditional expectation and variance are strictly linear functions of the state x. 158 The generalization to affine functions of x is provided in S3 Appendix. 159 Example Model Types 160 The following models satisfy the sufficient conditions: 161 Continuous-Time First-Order Reaction Networks: 162 A continuous-time first order reaction network is specified by a set of reactions R, 163 and a propensity λ r (x, t) and stoichiometry s r for each reaction r ∈ R. The 164 propensity λ r (x, t) is the expected rate at which reaction r occurs at time t. In a 165 first order reaction network, the propensities are linear functions of x, denoted p r x. 166 The actual reaction times are random. The probability that reaction r occurs once 167 in a small time interval is λ r (x, t)∆t + O(∆t 2 ), doesn't occur is 168 1 − λ r (x, t)∆t + O(∆t 2 ), and occurs more than once is O(∆t 2 ). If reaction r occurs 169 at time t then X(t) is replaced with X(t) + s r . The Langevin approximation to a continuous-time first order reaction network is 173 dX(t) = AX(t)dt + B(x, t)dW (t) where W (t) is a Wiener process, A = r∈R s r p r , 174 and B(x, t) = r∈R λ r (x, t)s r e r where e r is the r th |R| × 1 indicator 175 vector [7,12,13]. The associated diffusion matrix, D(x, t) = r∈R λ r (x, t)s r s r , 176 governs the instantaneous production of variance, while A governs the expected 177 behavior of X(t). with mean λ r (X t , t)∆t (c.f. [2]).

183
Discrete-Time Matrix Models: 184 A discrete-time matrix model is a model where the distribution of X t+∆t , given 185 X t = x, is parameterized by a linear function of x, say Ax for some matrix A, and 186 some fixed external parameters. Usually the expected value of X t+δt is given by Ax. 187 For example, we could set each entry of X t equal to a Poisson random variable, with 188 means fixed by AX t−∆t .

189
In practice, not all discrete-time matrix models use Poisson distributed random 190 variables. Instead, the product a ij x tj may represent the mean of a different distribution. In an age or stage-structured population model, the entries of the first 192 row of A represent fecundity rates (reproduction), and the remaining entries all 193 represent transition probabilities between stage/age classes j and i [14]. Then, the 194 actual number of individuals who transition from class j to class i is binomially 195 distributed, with x tj trials and probability a ij . Thus the expected number of 196 individuals making the transition is a ij x tj , but the variance in the number is 197 a ij (1 − a ij )x tj , not a ij x tj . Nevertheless, the variance remains a linear function of x. 198 Transitions representing survival, dispersal, aging, or growth often take this form.

199
Binomial survival distributions appear widely in discrete time population models 200 subject to demographic stochasticity [15].

201
Fecundity distributions are more model dependent, and may not be known 202 explicitly [16]. While some authors and software packages adopt the Poisson 203 distribution [17,18], the Poisson distribution is often ill-suited to modeling the 204 distribution of offspring since it often lacks biological justification, is unbounded 205 above, and may have unrealistic variance given its mean. While a different 206 parametric family of distributions could be used, for example, the generalized 207 Poisson recommended in [16], our analysis depends only on the functional form of

216
Theoretical Results

217
If the stochastic process satisfies conditions P1-P4 above, then the state covariance 218 takes a standard form. Let V (t) or V t denote the covariance in the state, andx(t) orx t 219 denote the expected state, at time t. Here we provide an explicit derivation for V t given 220 x 0 and V 0 . For simplicity, we assume that the time units are chosen so that ∆t = 1.

221
The results for continuous-time systems are entirely analogous. We articulate the 222 explicit forms for each case. Appendix S1 Appendix provides the derivation for 223 continuous-time systems. In this section we restrict attention to the case when the 224 conditional expectation and variance are strictly linear functions of x. The affine case is 225 treated in S3 Appendix, and results are reported at the end of that section.

226
Under conditions P1-P2, the expected state given the previous state is a 227 time-independent linear or affine function of the previous state. In the linear case there 228 exists some update matrix A such that: The covariance also satisfies a recursion. Since the process is Markovian, we can 231 apply the law of total variance to decompose the covariance at time t + 1: The expected value of X t+1 given X t is AX t+1 so: This term represents the forward propagation of past uncertainty into the present. The 234 product A(·)A is the forward propagation operator.

235
Under condition P3, there exists some linear function of x, D(x, t), such that: Then, since D(x, t) is linear: Therefore the covariance satisfies the recursive update equation: The recursion closes, leaving: The state covariance for a continuous-time process takes an analogous form: In continuous time the forward propagation of 248 uncertainty is governed by e At (·)e At , and uncertainty is produced continuously, so the 249 sum over the past is replaced with an integral.

250
In both cases, the second term is the variance produced by the process during the 251 time interval [0, t]. We will generally focus on this term, since the variance inherited 252 from uncertainty in the initial conditions is independent of the noise sources. The 253 variance produced by the process with initial expectationx(0) is the same as the 254 variance in the process if initialized at X(0) =x(0). Therefore, for the rest of this 255 discussion we drop the initial variance from Eq. 7 and Eq. 8.

256
Noise Source Expansion and Importance 257 Under condition P4, we can decompose the variance production term D(x, t) into a sum 258 over noise sources. Let N denote the set noise sources, and let D n (x, t) be the 259 variability in X t or X(t) produced by the n th noise source. Then: The variance produced by the process is: Lemma 2: (Additivity) If S(t) and S (t) are disjoint subsets of N at all t, then the state covariance in the process with noise sources S ∪ S equals the sum of the covariances in the reduced processes with noise sources S and S separately.
The second property of eq. 10 involves ordering of covariance matrices. Covariance 267 matrices can be (partially) ordered since they are real-symmetric and positive 268 semi-definite. If B and C are both Hermitian positive-definite then B ≥ C if B − C is 269 positive semi-definite, and B > C if B − C is positive definite [19]. If || · || is a 270 monotone matrix norm then B ≥ C implies ||B|| ≥ ||C||. From here on · will denote 271 a monotone matrix norm.

272
Eq. 10 ensures that the state covariance is monotonically nondecreasing as more 273 noise sources are added. The noise source expansion in Eq. 10 represents the variance 274 produced by the process as a sum of positive semi-definite matrices, one for each noise 275 source. Therefore:

276
Lemma 2: (Monotonicity) If S (t) ⊆ S(t) for all t and V (t) and V (t) are the covariances produced by the processes with noise sources T (t) and S (t) then V (t) ≤ V (t) and ||V (t)|| ≤ ||V (t)|| for any monotone matrix norm.
The importance of a set of noise sources can be analyzed by considering the variance 277 produced by a reduced process with the noise removed from the remaining sources. Let 278 S(t) denote the set of noise sources that are included in the process at time t. Then the 279 variance produced by the reduced process is given by restricting the sum over noise 280 sources to the sum over S.

281
To analyze noise source importance, we also need to specify how we measure the size 282 of a covariance matrix, and whether we are interested in the full state of the process or 283 some observable function of the state. Here we restrict our focus to observables that are 284 linear functions of the state. Let M be the Jacobian of the observable so that each 285 column of M corresponds to a specific observable, and inner products with columns of 286 M represent measurements. The matrix M is the measurement matrix. Then the 287 importance R of noise sources S at time t, given initial expectationx 0 , and 288 measurement matrix M , is the ratio of the norm of the covariance in the observable for 289 the reduced process, to the norm of the covariance in the observable for the full process. 290 That is, if V (S) t denotes the state covariance in the process with sources S, then: Then, for discrete and continuous-time respectively: Note that if the chosen norm is linear in the matrix entries, then the importance can 293 be interpreted as the fraction of the uncertainty contributed by noise sources S. It is 294 also natural to seek a matrix norm that is unitarily invariant, so that the size of the 295 covariance is independent of our representation of the state space. The only linear, 296 unitarily invariant matrix norm on symmetric-positive definite matrices is the trace 297 norm. The trace is a monotone norm for positive semi-definite matrices since the 298 diagonal entries of a positive semi-definite matrix are all nonnegative. The trace is also 299 a natural choice in this context, since the trace of a covariance is the expected (squared) 300 distance between a sampled observable and its expectation.   If A is diagonalizable then A = U ΛW where U and W are matrices whose columns 306 are the right and left eigenvectors respectively, and W = U −1 . Then: Next, since D n (x, t) is linear in x, and time independent: Therefore, the variance produced by the process at time t is: where the sum over noise sources can be moved outside the sum over time via the 310 assumption that the set of noise sources is unchanging.

311
Next, define the matrix: Then Eq. 15 can be expanded into a sum over triples of eigenvalues: The sum can be rewritten as a geometric series, so: In the special case when λ i λ j = λ k the sum equals λ k t t.

315
Therefore, the closed form for the covariance produced by the reduced process with 316 noise sources S is: The analogous equation in continuous-time is: Then, to compute finite time importance, substitute Eq 19 or 20 into Eq 11. Note 319 that the existing measure introduced by Schmidt and Thomas, only considered 320 importance in the long time limit [1]. Finite time importance is of interest in 321 applications whenever fluctuations about transients matter, such as during an invasion 322 process.

323
The only difference in the discrete and continuous-time expressions is the time  Eq. 11 establishes a closed form for finite-time noise source importance. The long-time 329 limit remains.

330
For processes with a non-degenerate steady state distribution, the state covariance 331 converges to a fixed matrix. For processes without a steady state, the covariance can 332 diverge or converge to zero. In these cases, we scale the covariance by its asymptotic 333 growth rate, and analyze importance with respect to the scaled covariance. Scaling the 334 variances by the asymptotic growth rate does not change the importance measure since 335 ||αB|| = |α| · ||B||, so the numerator and denominator are scaled by the same constant. 336 To analyze the long-time limit we require an additional property:

339
The limiting behavior of the variance depends on the limiting behavior of Eq. 18 as t 340 gets large. These limits depend, in turn, on whether |λ 1 | is greater than, equal to, or 341 less than one, which determines whether X t is expected to grow, remain constant, or 342 shrink. When growing, the variance grows proportionally with λ 2t 1 , when neutral the 343 variance grows proportionally with t, and when shrinking the variance shrinks 344 proportionally with λ t 1 .

345
The respective limits are:

346
Growing: if |λ 1 | > 1 then lim In each case, the scaled variances take the general form: where J is an index set and f (Λ) is a function of the eigenvalues. Both J and f (Λ) 348 depend on the limiting scenario, cf. Table 1. Table 1. Index set and eigenvalue dependence for long time limit of discrete-time variances, scaled by their growth rate.
Note that each of the cases has different limiting behavior. For example, when 350 growing, the covariance converges to a constant when rescaled by λ 2t 1 , the growth rate of 351 ||x(t)|| 2 . Then the limiting behavior of the scaled covariance is equivalent to the ||x(t)|| , the relative fluctuations about the expected state. When 353 shrinking, the variance is proportional to λ t 1 instead of λ 2t 1 . The difference in these 354 limiting behaviors can be explained by considering the primary sources of uncertainty at 355 long times.

356
If X(t) is expected to grow, then the expected trajectories diverge exponentially. As 357 a result, the uncertainty at some long time is mostly generated by uncertainty produced 358 at the start of the process. Thus, uncertainty early in the process (far in the past) is   What if X(t) is neither expected to grow nor shrink? Then there is a direction along 378 which X(t) is a martingale. Specifically, if λ 1 = 1, then the component of X(t) along 379 the dominant eigenvector, u 1 , is an unbiased random walk. This random walk accrues 380 variance linearly in time. Past variance does not decay or grow, so the variance along u 1 381 grows steadily with each added step. Meanwhile, the projection of X(t) onto all other 382 eigenvectors vanishes, so the long time variance is dominated by the linear growth along 383 u 1 .  Table 2.

387
The finite time noise source importance depends on the set of noise sources 388 considered, the measurement matrix (the choice of observable), the chosen norm, and 389 the initial expected state of the process. Our next Lemma answers the question: What 390 does the long time noise source depend on? Table 2. Index set and eigenvalue dependence for long time limit of continuous-time variances, scaled by their growth rate.
Lemma 3: In the growing case, the long time noise source importance, lim t→∞ R(S|M, || · ||,x 0 , t), is independent of the choice of measurement matrix M and norm || · ||. In the shrinking case, the long time noise source importance is independent of the initial conditionsx(0). In the neutral case without a conservation constraint, the long time noise source importance is independent of the choice of measurement matrix M , the norm || · ||, and the initial conditionsx(0).
Proof: In both the growing and neutral cases (without conservation constraint) all 392 contributions to the covariance take the form of a scalar coefficient times M u 1 u 1 M . It 393 follows that the variance associated with a set of noise sources is the sum of those scalar 394 coefficients times the matrix M u 1 u 1 M . The same holds when using all noise sources. 395 Therefore the covariance in any reduction of the process is proportional to the 396 covariance in the full process. The ratio of the norm of two proportional matrices is the 397 ratio of their proportionality constants, since ||αA|| = |α|||A|| for any norm. Thus, the 398 asymptotic importance in the growing and neutral case only depend on the scalar 399 coefficients that scale M u 1 u 1 M .

400
In particular, in the growing case, the asymptotic importance for discrete and 401 continuous models takes the form: In the shrinking and (unconstrained) neutral cases, all contributions to the 403 covariance are proportional to (w 1x (0)), so the dependence on initial conditions cancels. 404 Thus, in the shrinking case: and in the unconstrained neutral case: Therefore in the growing case the importance is independent of the measurement matrix 407 and choice of norm, in the shrinking case the importance is independent of the initial 408 conditions, and in the unconstrained neutral case the importance is independent of the 409 measurement matrix, norm, and initial conditions.

410
Lemma establishes that noise source importance is independent of what is observed 411 when the process grows or wanders. This result stands in contrast to the existing 412 literature, where the choice of measurement strongly influences which noise sources are 413 important. Indeed, the stochastic shielding heuristic as originally proposed by 414 Schmandt and Galan [2] is predicated on the assumption that the only important noise 415 sources are those that directly change the observable. In that case, edge importance is 416 entirely determined by what is measured. Schmidt and Thomas introduced their 417 measure to show that the heuristic does not always hold [1]. Nevertheless, edge 418 importance reversal usually occurs when the dynamics of the observable are strongly 419 coupled, albeit indirectly, to a distant noise source [7]. Thus, importance is still directed 420 by the observable.

421
As Lemma ?? establishes, the long time noise source importance does not depend on 422 the choice of observable, when the triple sum over i, j, k is restricted to a sum over k.  17)). This occurs when most of the variance is produced in 426 the distant past, as in the growing case, or when it accrues linearly along one direction 427 and shrinks along all others, as in the wandering case. Note that, in both cases, the produces variance along w 1 , so the inner product of w 1 X t is conserved.

446
Thend 1,: (k, n) = 0 andd :,1 (k, n) = 0, since w 1 is orthogonal to the range of D, so 447 any term in the triple sum in Eq. 17 with i = 1 or j = 1 is automatically zero, which 448 sets the diverging term to zero. If k = 1 then the entire term converges to zero, so the 449 only surviving terms are k = 1, i = 1, j = 1. Then the long time limit (steady-state) 450 covariance equals: where C = w 1x 0 = w 1 X(t) is the conserved quantity.

452
The analogous result in continuous-time is: To compare to the edge importance measure proposed in [1] we write outd i,j (1, n) 454 in detail. First, note thatD(k, n) = W D n (u k )W. Next, if the continuous-time model is 455 a reaction network, then each reaction r is a noise source, so we replace n with r and 456 D r (u 1 ) = (p r u 1 )s r s r . Then: In a continuous-time first-order reaction network u 1 is proportional to the 458 steady-state distribution, and Cp r u 1 is the rate at which reaction r occurs at steady 459 state. This flux, J r , is the steady state flux of reaction r. Then the uncertainty 460 contributed by a single reaction is: The associated asymptotic importance is: which coincides with the edge-importance measure defined by Schmidt and Thomas for 463 first-order continuous-time reaction networks when the observable is one-dimensional 464 (M has only one column) [1].

466
The previous sections assumed that the conditional moment equations were linear 467 functions of the state. Here we generalize our results to the affine case. The derivation 468 largely follows from the linear case and is provided in S3 Appendix.

469
In the affine case Eq. 1 and Eq. 4 are replaced with: where each D n (x, t) is an affine function of the state x, D  Alternatively, if the linearization is based on a best fit linear model to observed 490 dynamics then the linearization will include state-independent terms whenever the fit 491 intercepts are nonzero. We consider an example of this kind in §Wood Frogs.

492
If the conditional moment equations are affine, then the moments satisfy the same 493 recursions as before,x t+1 = Ax t + b t , and V t+1 = AV t A + n∈N D n (x t , t). The 494 recursions are unchanged since they depend only on the linearity of expectation, law of 495 total variance, and linearity of the conditional moments.

496
While the recursions are unchanged, the long time limits of the moments are 497 changed by adding state independent terms. In particular, adding constant source terms 498 to the conditional moment equations introduces a new limiting scenario.

499
In absence of a constant source term, the only limiting scenario with a 500 non-degenerate steady state distribution was the neutral case, λ 1 = 1, with a 501 conservation constraint. In the shrinking case, |λ 1 | < 1 the steady state distribution was 502 a delta distribution at the origin. In the presence of constant source terms the |λ 1 | < 1 503 case admits non-degenerate steady states. If b > 0 and |λ 1 | < 1, then the expected state 504 will converge towards a nonzero equilibrium, and thus the variance source terms will 505 approach constant, nonzero values. The resulting steady state distribution balances the 506 fluctuations produced by the nonzero noise sources, with their tendency to decay back 507 towards the equilibrium. This is the type of steady state commonly observed in 508 Ornstein-Uhlenbeck processes. If b = 0 thenx t will converge to the origin. In the linear 509 case the intensity of each noise source is proportional to x, so whenx t approaches the 510 origin all the noise sources stop producing noise -hence the approach to a delta 511 distribution. In the affine case, when b = 0 it must be true that D Wandering: if λ 1 = 1 then: 1,1 (k, n) + . . .
HereD = W DW ,x denotes expansion on the eigenbasis,x = W x, and x * denotes 520 the equilibrium of the expected dynamics, lim t→∞ E[ Once the long time variance has been computed, the importance of each noise source 522 can be analyzed using the usual approach. Computing the long time variance in the 523 affine case is considerably more expensive than in the linear case, as all scenarios 524 require computing a sum over each eigenvalue, and the steady state scenario requires 525 computing a sum over all triples of eigenvalues.

527
The four models introduced in §Scope differ only in the formulation of A and D n (x).

528
For both continuous-time reaction networks A = r∈R s r p r and D r (x) = (p r x)s r s r .

529
For a discrete-time reaction network with time step τ the matrices are 530 A = I + τ r∈R s r p r and D r (x) = τ (p r x)s r s r . For a discrete-time matrix model, A is 531 explicitly fixed by the model parameterization, while D n (x, t) will depend on the 532 distributions used to sample the updates. (100-700 µM) [28]. Calcium is released from the ER by second messengers to raise 554 intracellular Ca 2+ concentrations and trigger a cellular response. Inositol trisphosphate, 555 InsP 3 , is a ubiquitous second messenger involved the calcium-mediated processes listed 556 above. It triggers Ca 2+ release by binding and activating ion channels in the ER 557 membrane [22,29]. Sensitivity to InsP 3 is modulated by the concentration of free 558 intracellular calcium through calcium-induced calcium-release (CICR). CICR amplifies 559 the signal so that a few InsP 3 binding events can trigger a rapid release of calcium from 560 the ER into the cytosol. Free calcium then diffuses through the cytosol. If enough 561 calcium is released, then the diffusing calcium binds to an activating channel receptor 562 and triggers the release of more calcium producing propagating waves [29]. Calcium also 563 binds to an inhibitory channel receptor which inhibits further channel openings when 564 the free calcium concentration is high [30]. Then the free calcium is gradually 565 reabsorbed into the ER.

566
The combination of high InsP 3 sensitivity at low Ca +2 concentrations, inhibition at 567 high concentrations, and diffusion, allows complex temporal and spatial dynamics.

568
These dynamics range from stochastic "blips" released by individual InsP 3 channels, to 569 "puffs" released by a cluster of channels, to travelling waves [28]. Shuai  One of the nine subunit states comprises an "active" configuration. The full InsP 3 576 channel model contains four identical subunits, and the channel opens if at least three 577 of the subunits are in the active state. As shown in Figure 1, each subunit has three 578 binding sites, one for InsP 3 , and two for Ca +2 . One calcium binding site activates the 579 channel, and the other inhibits the channel. The first is responsible for CICR, while the 580 second is responsible for closing the channel when free Ca +2 concentrations are high.

581
Following [30], we use a 1 to denote a bound site and a 0 to denote an unbound site. probability that the channel is active, and the expected time the channel spends closed 590 as a function of calcium concentration [30]. The exact model and parameters are 591 provided in [30]. Note that all binding events are first order reactions, and occur with 592 rates that depend on the InsP 3 and Ca +2 concentrations.
Since the model conserves the total number of subunits, it admits a steady state 596 distribution and falls into the fourth limiting scenario ("persisting"). We study the 597 importance of each noise source to the steady state variance in the number of open 598 channels as a function of the InsP 3 and Ca +2 concentrations. All importance values are 599 calculated using Eq. 29. Finite time importance can also be computed using Eq. 11.

600
Finite time importance could be relevant when considering the channel's response to a 601 spike in InsP 3 or Ca +2 . We focus on steady state importance here.

602
To compute importance, we specify a set of noise sources. We start by considering 603 individual edges as noise sources; however the edges could also be grouped by reaction 604 type. The importance of a reaction type can be recovered by summing the importance 605 of the edges in that class. Four pairs of binding and unbinding edges correspond to the 606 InsP 3 binding site. Similarly, four pairs of edges correspond to binding and unbinding of 607 the activator calcium binding site, and four pairs correspond to binding and unbinding 608 at the inhibitory calcium site. These sets correspond to all edges of the cube that move 609 up and down (InsP 3 ), forward and back (activator calcium), and left and right  Importance is not evenly distributed among the 13 reaction pairs. For any 627 combination of InsP 3 and Ca +2 concentrations, the 4 most important edges account for 628 more than half the steady state variance, and the 8 most important account for more 629 than 90% of the steady state variance. The number of reactions needed to account for 630 50%, 90%, and 95% of the steady state variance is shown in Figure 2. Three of the four 631 edge pairs involving the active calcium binding site never contribute more than a percent 632 of the total variance. Binding of activating Ca +2 only matters if the inhibitory site is 633 unbound, and the InsP 3 site is bound. Most of the edges contribute less than a percent 634 of the total variance when the calcium concentration is low (< 0.1 µM) or the InsP 3 635 concentration is low. Consequently, most edge pairs do not contribute significantly to 636 the steady state variance unless the calcium concentration is intermediate (1 to 10 µM) 637 and the InsP 3 concentration is small (< 1 µM). Similar results extend to the classes of 638 edges corresponding to particular reactions. Figure 3A shows regions in the Ca +2 InsP 3 639 plane where each class of reactions is the least important.

640
Since most of the reactions do not contribute significantly to the variance in the 641 number of active subunits the simulation of the InsP 3 channel could be streamlined by 642 simulating unimportant reactions deterministically rather than stochastically. If the free 643 Ca +2 and InsP 3 concentrations are held fixed for the duration of the simulation, then 644 the importance of each reaction could be computed a priori, and unimportant reactions 645 could be identified and simulated deterministically. This approach is easily implemented 646 for Langevin approximations (c.f. [1,7]) to the reaction network, and reduces the Fig 2. The number of edges needed to reach 50%, 90%, and 95% of the total variation in the number of active subunits depends on the InsP 3 concentration (marked I) and the free Ca +2 concentration (marked C). All concentrations are in µ M and range from 10 −2 to 10 1.75 µ M. We never need more than 19 of the 26 edges to reach 95% of the total variance, and do not need most of the edges unless the Ca +2 concentration is intermediate and the InsP 3 concentration is low. computation cost by reducing the total number of random numbers needed to perform 648 the simulation.   Figure 3A. In contrast, the importance of the observable 673 reactions does not approach zero anywhere, and is never less than 0.14 (high Ca +2 , low 674 InsP 3 concentrations). The observable reactions are never of negligible importance since 675 the observable only changes when observable reactions occur. Thus, even if other 676 reactions are more important, the observable events still matter.  The linearization produces affine moment update equations, so requires the fully general 692 analysis presented in §Affine Conditional Moments. We consider importance at both 693 long times and at short times during a recolonization process. The bottom row shows each surface as a heat map. All panels share the same color bar. The red lines mark importance 0.5. Note that the importance of activation is largest when the Ca +2 concentration is small, and decreases as the concentration increases. Binding to the active site goes beneath 0.5 at the boundary where the green region ends in Figure  3B. Binding to the active site becomes unimportant when the Ca +2 concentration is large, near the boundary of the yellow region in Figure 3B. The importance of the inhibitory reaction is large when the Ca +2 concentration is large, and decreases as it decreases. Like the activation reaction, the inhibition reaction's importance falls beneath 0.5 at the boundary of the yellow region, and becomes negligible in the green region of Figure 3B. This leaves an intermediate Ca +2 concentration regime (between 0.1 and 10 µM) where the other two reaction classes are important. The InsP 3 binding reaction is important when the InsP 3 concentration is low, while the observable reactions (conformational change in the subunit between active and inactive) is important when the InsP 3 concentration is large. The importance of the observable reactions closely follows the probability that the channel is active. of distributions. Where needed, we estimate conditional variances using data in [33]. 706 Krüger's model is an age-structured model containing 37 age classes (ages 0 through 707 36 years). All individuals of age 36 are assumed to die. Eagles begin breeding after age 708 four. Separate life history rates were estimated for the periods 1947 to 1974, when the 709 species was endangered due to pollution, and 1974 to 2008, when the species recovered. 710 We assume that eagle deaths are independent, and that within a given time period 711 (1947-1974 or 1974-2008), the probability an eagle dies depends only on its age. Then 712 the number of eagles who survive in a given age class each year is a binomial random 713 variable, with parameters fixed by the survival rate for that age class, and the number 714 of individuals in that age class. The expected number of surviving birds from each age 715 class is linear in the number of birds in the preceding age class, as is the variance. This 716 model is supported by [34] who found no evidence for density dependence in eagle 717 survival rates.

718
Survival rates were estimated by [33] based on age of death for 80 birds between 1953 719 and 2008. Due to small population sizes only 25 dead birds were recorded before 1974. 720 As a result, survival rates were only estimated up to age 4 for the 1947-1974 period, and 721 all survival rates for older eagles were set equal to the estimated survival rates for 722 1975-2008 [33]. The impacts of DDT on egg, chick, and adult survival rates are studied 723 in [32,35,36]. Not all survival rates could be estimated directly. Survival rates for age 724 classes missing data were set to the average of survival rates for neighboring age classes. 725 Here we present results based on Krüger's reported survival parameters [33]. These 726 parameters appear noisy, likely due to limited sample size. In order to evaluate the 727 impact of parametric noise on our results, we also considered a version of the model in 728 which survival rates as a function of age are smoothed after age 5. Survival rates were 729 only smoothed after age 5, for two reasons. First, chick mortality is higher than juvenile 730 mortality. Second, mortality at age 4 is high due to territorial competition between 731 immigrant eagles of age 4 seeking a nesting site and adult eagles with established nests 732 who defend their territories [34,37]. We found that the smoothed model generated 733 smoother trends in importance across ages, but as our qualitative conclusions were 734 unaffected, we present results for the original unsmoothed model here. 735 We modeled the birth process as follows. The probability a particular female of a 736 given age breeds is provided by [33]. Breeding probability is lower for young adults 737 (0.620 at age 5) and approaches 1 for older adults (0.970 by age 8). Assuming  To study the importance of each noise source we need an estimate for the variance in 748 fecundity (fledged chicks per breeding female). Unfortunately, this variance is not 749 provided explicitly in [33]. We used a two-pronged approach to bridge this gap. On one 750 hand, we used clutch size data from [33] to estimate the variance. We validated our 751 estimates against studies of Lithuanian and Swedish white-tailed eagle populations 752 (see [38,39]). On the other hand, we developed a variational approach that provided 753 interval estimates for the importances. In the variational approach, we impose 754 biologically motivated constraints on an unknown distribution, fix its mean, then 755 maximize and minimize its variance given the constraints. We use this approach to 756 avoid introducing specific distributional assumptions that are convenient but lack good 757 biological motivation. For example, we avoid assuming that clutch sizes are Poisson 758 distributed, as is assumed in some models (c.f. [34]), since clutch sizes above three or 759 four eggs are virtually never observed and are physiologically implausible. Instead, we 760 assumed that eagles: a) lay at most four eggs, b) lay zero eggs with probability 0.7 for 761 1947-1974 and 0.2 for 1975-2008 [33], and c) lay 3 eggs more frequently than they lay 4 762 eggs, as was observed in empirical clutch size distributions [38,39]. Then, for each age 763 class, we found the distribution with largest and smallest possible variance satisfying 764 constraints a-c with mean equal to the average fecundity provided in [33]. The 765 variances estimated from clutch size data [33] always fell between the max and min 766 variance provided by the variational approach. Therefore, we report results from the 767 variational approach which demonstrate the range of plausible importances across the 768 range of plausible fecundity variance.

769
Eagle sex ratios at birth are close to 0.5 [33], so we assume that the number of 770 female fledged chicks is a binomial random variable with mean equal to half the number 771 of fledged chicks. We do not model males.

772
Each step in the life cycle is a separate noise source. We can then decompose the 773 variance produced by birth into components associated with separate life events. Since 774 we only track the number of fledged female chicks produced at the end of the birth 775 sequence, we use the law of total variance to decompose the variance in the number of 776 fledged female chicks into contributions from each random event (whether a female 777 breeds, clutch size, and sex determination). The sum of these contributions is the Krüger's predictions [33].

793
To compute noise source importance, we set the initial population equal to ten 794 individuals, and partition the population into age classes according to the stable age 795 distribution. We then compute the importance of each noise source to the total 796 population using equations 23 and 24. Figure 5 illustrates results for the maximum 797 variance case.

798
In all cases, we find that mortality contributes more to the overall noise in eagle 799 populations than reproduction, with mortality accounting for 59 to 62% of variance 800 before 1974 (max and min variance in clutch sizes), and 79 to 82% of all variance after 801 1974. After 1974 the importance of mortality peaked at age 4. This peak can be 802 explained in two ways. First, eagle mortality is high at age 4 due to territorial contests 803 between young adults hoping to establish a new nest, and established adults. Notably, 804 this increase in mortality is not observed before 1974 when eagle populations were small, 805 so nesting space was less limited [33]. Second, the transition from age 4 to age 5 is 806 crucial, as it represents the transition to sexual maturity. Before 1974 the most 807 important noise source is chick mortality. We find that, across most population models 808 studied, the most important noise source is either mortality preceding the transition 809 into sexual maturity, or mortality during the first year of life, which is often high. In 810 both eras, we see large noise source importance for mortality during the first 8 to 10 811 years of life. These are crucial years, since individuals must survive the first four years 812 to reproduce, have notably lower survival rates as young adults, and take approximately 813 8 years before they breed reliably.

814
Reproduction accounts for 38 to 40% of the total variance in population size before 815 1974 (min and max variance in clutch sizes), and 18 to 20% of the total variance after 816 1974. The variance contributed by whether or not a female breeds is largest at age 5 817 and declines quickly as 97% of all females older than 8 breed. Thus, uncertainty in 818 breeding generally contributes the least variance of all noise sources considered (≤ 1%), 819 and is negligible for eagles over 8 years old. Before 1974, variance due to clutch size is 820 consistently larger than variance due to individuals' sex, accounting for 26 to 29% of the 821 total variance. Variance due to individuals' sex accounted for 11 to 12% of the variance 822 in these cases. After 1974, clutch size and sex contribute comparable amounts of Top row: Importance of noise sources for variance in white-tailed eagle populations before 1974, assuming maximum variance in clutch sizes. Note the high variance from chick mortality, large noise contributions from ages 4 and 8, and the slow decay in the importance of reproduction. The erratic behavior of importance in mortality at later ages is likely due to noise in the estimated survival rates, and was smoothed out after smoothing the survival rates (not shown). Bottom row: Importance of noise sources for white-tailed eagles after 1974, assuming maximum variance in clutch sizes. Note the spike in importance at age 4, preceding the transition to sexual maturity, and for young adults experiencing greater mortality.
of each reproduction event peaks at age 6, and then declines steadily with increasing age. 825 Accounting for all sources of noise at a given age we find that, before 1974, noise Wood frogs are a well-studied (c.f. [41][42][43]) North American amphibian. In order to 855 study the relative importance of distinct noise sources within wood frog populations, we 856 linearize a nonlinear metapopulation model. We consider two distinct ecological 857 scenarios: equilibrium across life stages and local habitats, and global recovery from a 858 large perturbation. Because the linearized conditional moment equations are affine, we 859 apply the results from §Affine Conditional Moments. In addition, this example 860 demonstrates how to analyze noise source importance along transients, in particular 861 during a colonization process following a large perturbation.

862
Adult wood frogs migrate to ponds to breed in the early spring, where adult females 863 lay a single egg mass. There, the eggs hatch into fully aquatic larvae which 864 metamorphose into terrestrial juvenile frogs in the early summer. Wood frogs typically 865 reach sexual maturity after two years on land before they return to the pond to breed as 866 what we will refer to as "young adults". Some frogs, which we call "mature adults", 867 survive to breed a second time in the following year. When wood frog ponds are spread 868 across a large terrestrial habitat, they may act as a metapopulation [42,44]. While the 869 majority of individuals return to breed in the pond they hatched in, some frogs disperse 870 during their juvenile life stage and breed in a new pond [42]. Variation in dispersal egg to juvenile is density-dependent (due to larval competition within each pond), as is 882 survival from juvenile to young adult (reflecting juvenile competition across the whole 883 terrestrial environment) [43,45,46]. Each of these density-dependent transitions uses a 884 Beverton-Holt model [47], which outperformed linear versions of each transition in 885 AIC [48] goodness of fit tests. Note that the model ignores certain important features of 886 the ESGR by assuming that the ponds are equivalent, thereby ignoring the differences 887 in predation between large and small ponds. For these reasons, our results should not 888 be interpreted as specific predictions about the ESGR, which would require a more 889 detailed model, but as illustration of our theory in a model with realistic age, stage, and 890 spatial structure.

891
To make the model amenable to our measures, we replace it with a stochastic model 892 whose conditional moment equations are affine. We perform this conversion in two steps. 893 First, we define a fully stochastic model consistent with the deterministic nonlinear The stochastic model follows. First, all young and adult female frogs reproduce in 900 each pond. While the form of the egg mass distribution is not known (roughly 901 log-normal), the necessary moments can be estimated from available data. Each female 902 lays 723.7 ± 197.3 eggs (mean and standard deviaiton). For the sake of simulation, we 903 use a log-normal distribution, rounded to integer values. We assume that the number of 904 eggs laid by each female is independent, so the variance in the total number of eggs laid 905 in a pond is the sum of the variance in the number of eggs laid by each female. We do 906 not consider sex at birth, or the probability a female does not reproduce, because data 907 on tadpole sex is not available and females who did not return to the pond habitat to 908 breed cannot and need not be distinguished from females who had died. Egg survival to 909 the juvenile stage is binomial with a density dependent survival probability: where s 1 is the survival probability of a single egg in the absence of any other eggs, M proportional to exp(−d/d 0 ) where d 0 is a reference distance (average dispersal distance 921 of frogs that do not return to their natal pond). Note that this dispersal model does not 922 distribute frogs equally across the ponds. Instead, well connected ponds tend to receive 923 more immigrants than they produce emigrants, so have higher total steady state 924 populations. Since survival is nonlinear and the ponds don't maintain the same total 925 populations, the ponds also don't share exactly the same age-distribution, though 926 dispersal is rare enough that the effect on the age distribution is small.

927
Once a breeding pond has been assigned to each juvenile, we sample a fraction that 928 survive to arrive in that pond at the end of the two-year juvenile stage. Survival 929 probability s L is density dependent, and depends on the total number of individuals in 930 the terrestrial environment. Specifically, where L is the total lag population (individuals in the terrestrial environment that have 932 not yet reached breeding age), s 1 is the survival probability for an individual in absence 933 of any others, and M 2 is a half-saturation constant. The number of young adult 934 individuals that survive to breed in pond i is then sampled from a binomial distribution 935 with survival probability set to s L (L) and number of trials set to the number of lag 936 individuals assigned to migrate into pond i.

937
Young adults who survive another year return to their breeding pond one year later 938 as mature adults. Survival of the young adults is not density dependent, so the number 939 of mature adults each year is drawn from a binomial distribution with a fixed survival 940 probability s Y . All mature adults die in the subsequent year.

941
To validate the stochastic model, we ran the deterministic ODE model to 942 equilibrium, and compared the equilibrium to the average populations in each age class 943 and pond after a long (2 × 10 5 years) trial run of the stochastic model. We found that 944 the deterministic equilibrium agrees with the average predicted by the stochastic model, 945 to within a maximum relative error 0.6% over all age classes and ponds.

946
The stochastic model does not have linear conditional moments, because egg and lag 947 survival are both density dependent (see equations 33 and 34). So, to apply the noise 948 source importance measures, we require a linearization procedure. We linearize by   The number of young adults who arrive at pond k depends on both the survival rate 959 of frogs during the pre-reproductive ("lag") terrestrial phase, which depends on the total 960 lag population, and the number of surviving lag individuals who choose pond k. Thus, 961 we fit the expectation and variance in the number of immigrating lag individuals into 962 each pond to the number assigned to that pond and the total number of lag individuals. 963 January 6, 2022 27/35 The associated fits were E[Y t+1 (j)|L t = l] = 0.031l(k) − 0.00065 i l(i) +  First, we observe that the importances depend much more strongly on which life 980 event is considered, than which pond it occurs in. Figure 7 illustrates this effect: the 981 variance in importance between life history events (differences among box heights) is 982 much larger than the variance across ponds (size of box and whisker associated with 983 each event).

984
Broadly speaking, the most important event is juvenile survival, followed by young 985 adult survival, egg survival, egg count, then dispersal, though the importance order may 986 change depending on which life stages are observed. For example, egg survival and 987 count are more important to variance in juvenile population size than young adult 988 survival is, but young adult survival is more important than any other event, even 989 juvenile survival, to variance in the mature adult population. Juvenile survival is the 990 most important event in most cases since it is the survival stage preceding reproduction. 991 Dispersal is the least important in all cases since, in the linearized model, the ponds are 992 indistinguishable and interact only with the total lag population, which is independent 993 of pond assignment. Consequently, while dispersal contributes to the variance in the 994 number of frogs in any individual pond, it does not contribute any variance to the total 995 population across ponds. Dispersal also plays a small role since the 94% of all frogs 996 return to their birth pond.

997
So, dispersal (and the associated spatial component of the model) is not an 998 important source of randomness in the process if we consider the total population across 999 all ponds after reaching a stationary distribution. However, if we consider specific 1000 ponds, or a transient away from equilibrium, then dispersal plays a more important role. 1001 For example, if only a single pond is considered, then all five noise sources associated 1002 with that pond are more important than noise generated in other ponds, and match the 1003 ordering of importance shown for the total population across ponds. The importance of 1004 noise sources in the other ponds is generally very small. If we track the population in 1005 pond i, then the most important event in pond j is always dispersal, since the only way 1006 noise can propagate from pond j to pond i is via dispersal. This fact is an elegant 1007 illustration, in an ecological context, of stochastic shielding [1,2,7].

1008
Dispersal is also important during transients. During particularly dry summers, 1009 some of the ponds in the reserve may desiccate, killing off the aquatic life stages. In 1010 subsequent, wetter years, frogs will disperse out from the larger ponds to repopulate 1011 smaller ponds that dried out. To complete our study, we considered the importance of 1012 each different noise source over the course of an recolonization from the largest pond in 1013 the reserve. The largest pond in the reserve is Fish Hook Marsh (Pond 8 in Figure 6). 1014 We initialize the system by setting the frog population to 0 in all ponds except pond 8, 1015 and set the population in pond 8 to its multi-pond equilibrium value. This initial 1016 condition is far from equilibrium, so we need to relinearize the model. S4 Appendix 1017 describes the relinearization technique.

1018
After linearizing about the perturbed initial conditions, we computed the importance 1019 of each noise source to the total population in a set of target ponds. Figure 8 shows the 1020 results when the source pond is set to pond 8 and the target pond is set to pond 3. This 1021 example is representative. Pond 3 was chosen as the displayed target since it is 1022 connected to Pond 8 by a chain of smaller ponds, so, for the right dispersal parameters 1023 we might expect to see a cascade of importance in dispersal down the chain from 8 to 3. 1024 Pond 3 is also an interesting target since it is large, but far removed from 8, and 1025 neighbors a cluster of smaller ponds (16,5,18). 1026 In the first year, the most important noise source to variance in the target pond's 1027 population size (pond 3) is dispersal out of the source pond (pond 8). Dispersal out of 1028 the source pond remains among the most important noise sources in later years, but 1029 decays quickly in importance after the first year. Juvenile survival in the target pond 1030 replaces dispersal out of the source pond as the most important source in year two, then 1031 remains the dominant source of noise in all but the first year. As the frog population in 1032 the target pond grows, the importances of the other life events in the target pond grow, 1033 approaching the importance order seen at equilibrium. While the noise sources in the 1034 target pond are the most important sources at relatively long times, dispersal out of the 1035 source pond remains more important than most of the noise sources in the target pond 1036 (excluding juvenile survival) until 15 -18 years after invasion. This slow decay in the 1037 importance of dispersal is noteworthy given how little dispersal, contributes at 1038 equilibrium.

1039
The slow decay of the importance of dispersal out of the source after the first year is 1040 a result of the cumulative nature of the measure. Dispersal importance out of the source 1041 pond includes all variance contributed by dispersal out of the source pond up to the current year. Thus, since dispersal out of the source pond is necessary to seed the frog 1043 population in all other ponds, dispersal out of the source in year one is a crucial noise 1044 source that contributes a significant portion of the noise in the target pond at future 1045 times.

1046
Juvenile survival in the source pond is also important for similar reasons. Dispersal 1047 is a relatively rare event, so the target pond may not be seeded during the first year 1048 post invasion. In that case, survival of the juveniles in the source pond is important for 1049 seeding the target pond in subsequent years.

1050
The remaining important noise sources are all associated with the cluster of ponds 1051 neighboring the target pond (16,5,18). Juvenile survival in ponds 16, 5, and 18 are the 1052 next most important noise sources, and are ordered in importance by their distance 1053 from the target pond. Thus, spatial structure influences importance during the initial 1054 phase of colonization. This effect stands in contrast to the situation at equilibrium, in 1055 which spatial structure plays no role.
model with a linear model. The methods are compared, and their respective advantages 1109 and disadvantages are discussed. 1110