Smoothing in linear multicompartment biological processes subject to stochastic input

Many physical and biological systems rely on the progression of material through multiple independent stages. In viral replication, for example, virions enter a cell to undergo a complex process comprising several disparate stages before the eventual accumulation and release of replicated virions. While such systems may have some control over the internal dynamics that make up this progression, a challenge for many is to regulate behaviour under what are often highly variable external environments acting as system inputs. In this work, we study a simple analogue of this problem through a linear multicompartment model subject to a stochastic input in the form of a mean-reverting Ornstein-Uhlenbeck process, a type of Gaussian process. By expressing the system as a multidimensional Gaussian process, we derive several closed-form analytical results relating to the covariances and autocorrelations of the system, quantifying the smoothing effect discrete compartments afford multicompartment systems. Semi-analytical results demonstrate that feedback and feedforward loops can enhance system robustness, and simulation results probe the intractable problem of the first passage time distribution, which has specific relevance to eventual cell lysis in the viral replication cycle. Finally, we demonstrate that the smoothing seen in the process is a consequence of the discreteness of the system, and does not manifest in system with continuous transport. While we make progress through analysis of a simple linear problem, many of our insights are applicable more generally, and our work enables future analysis into multicompartment processes subject to stochastic inputs.


Introduction
Many biological processes comprise multiple independent stages.Viral replication, for example, is a multistage process whereby virions enter a cell through endocytosis, are unpackaged before DNA replication, repackaging, and release (fig.1a) [1][2][3].Similar multistage processes are evident in bacteriophage replication [4] and progression through the cell cycle [5], pervasive at the molecular (i.e., cascade reactions [6]) and macroscopic (i.e., transport through discrete layers [7]) levels, and even manifest in social processes such as queuing [8].A challenge for many systems is to modulate the impact of what are often highly variable external environments.For instance, while the intermediate stages of viral replication may be optimised to achieve high levels of virion multiplication, the system has either no, or only very limited, control over the number of virions entering the cell [9,10].For lytic viruses, should the number of virions present inside a cell exceed capacity the cell will lyse, destroying the system and ceasing replication [11].
The time until cell lysis-more broadly, the time until the first occurrence of any event within a stochastic process-can be modelled as a first passage time (FPT) and is dependent, among many other factors, on the variability and the autocorrelation of the process [12,13].
Statistics such as the variance, autocorrelation function, and the FPT, are commonly studied in scalar stochastic systems in biology [12,14].Many linear and non-linear systems described by continuous or discrete-space random walks have closed form solutions available for the aforementioned statistics, and if not for the FPT distribution itself, then for the mean, variance, and higher order moments of the FPT [12,15,16].Analysis of higher-dimensional systems (i.e., described by more than one first-order differential equation) has, to date, been restricted to a fixed number of dimensions; most commonly two or three, where the velocity or acceleration of a particle is described by a stochastic process [17,18].General techniques for analysis, such  We study a linear analogue of the virus problem, namely a system subject to random input, I(t), modelled as an Ornstein-Uhlenbeck process with mean µ, reversion strength θ, and noise magnitude σ.A realisation of the input is shown in (c).Material then progresses through n compartments at constant rate k.(d) A realisation of the system in (b), where X ν (t) models the concentration of material in each compartment ν = 1, 2, . . ., 6.It is evident that passage through the compartments has a smoothing effect.
as through the Fokker-Planck equation, quickly suffer from the curse of dimensionality, and dimension reduction techniques may yield lower-dimensional stochastic processes that lack the Markov property typically exploited in analysis.
Motivated in particular by the viral replication cycle, in this work we study the properties of a linear n-compartment model subject to an independent external input, which we model using a mean-reverting Ornstein-Uhlenbeck process (fig. 1b).While the presentation and analysis of stochastic models is ubiquitous throughout the biological literature, efforts have been largely restricted to the study of the effects of intrinsic noise, process noise, and fluctuating parameter values, including for analysis of viral replication dynamics [19][20][21][22][23].Meanwhile, there has only been a very small amount of work focussed on analysis of the statistical properties of largely deterministic systems subject to independent stochastic input [24].In the context of control and state-space estimation, our model parallels highly studied linear state-space models such as autoregressive models and the linear Kalman filter [25].Given the scarcity of general knowledge related to the behaviour of multicompartment models subject to input noise, we are motivated to study a simple linear analogue of the viral multicompartment problem subject to a single Ornstein-Uhlenbeck input.
The choice to study a simplified linear stochastic model allows us to formulate the multicompartment problem as a multidimensional Gaussian process, enabling us to draw on the significant body of literature devoted to study of the statistical properties of such systems [26][27][28] to formulate a series of analytical expressions for key statistics unique to stochastic processes including the variance, covariance, and autocorrelation function.Alternative approaches to a more theoretical analysis could include the study of system response to pure waves and input pulses; however, the goal of this work is to study the simple model directly.While we are not able to solve explicitly for the probability density function of the FPT, we present a series of numerical and approximate results that provide insight into the FPT, and the rate at which the mean FPT scales with the number of compartments in the system.We then apply our linear model to study how the behaviour or robustness biological systems can be modified through perturbations to unidirectional progression through the system.Viral replication, for example, is known to be a highly stochastic process, and progression through replication stages is very often not unidirectional [2].

Mathematical model
The problem presented in fig.1b can be expressed as the linear system of stochastic and ordinary differential equations Here, we denote by I(t) the random input, modelled as an Ornstein-Uhlenbeck process with mean µ, noise magnitude σ, and reversion strength θ; by X ν (t) the concentration of material in compartment ν; and by k the progression rate of material from one compartment to the next.All parameters are real and positive, and W (t) represents a Wiener process such that W (t + η) − W (t) is normally distributed with mean zero and variance η.
More compactly, we write where we define for X(t) = I(t), X 1 (t), . . ., X n (t) ⊺ .For notational convenience, we interchangably refer to I(t) as X 0 (t) (i.e., the input is thought of as compartment ν = 0).Equation ( 2) demonstrates that the system is a multidimensional Ornstein-Uhlenbeck process and, therefore, a Gaussian process.
We refer to Θ as the connectivity matrix, as it plays a role similar to that in graph and network theory, defining the connectivity between compartments in the system.Therefore, provided the system remains linear, we can arbitrarily express systems with any network structure (i.e., with non-local feedbacks or multiply-connected components) using eq.( 2).The form of µ, which contains the lower block matrix Θ 22 corresponding to Θ with the first row and column excluded, simplifies for the system in fig. 1 and eq.( 1) to µ, µ/k, . . ., µ/k ⊺ .Unless otherwise stated, we fix θ = µ = k = 1 and σ = 0.5 as default parameter values when producing simulation results.
There are various initial conditions that we consider in this work, based on the assumption that a stationary limiting distribution for X(t) exists (since θ, k > 0, this is always true for the form of Θ expressed above, and more generally provided that all eigenvalues of Θ have positive real part [26]).The first relevant initial condition is where X 0 is entirely specified.
We refer to this choice as the fixed initial condition.For the virus-cell lysis problem, we may be interested in setting X 0 = µ, 0, . . ., 0 ⊺ ; i.e., the concentration is zero in all compartments, and the input is initiated at its mean.A second, more biologically realistic initial condition, is where all compartments are initiated with zero concentration, but where the input is initiated from its stationary distribution I(0) ∼ N (µ, σ 2 /(2θ)).We refer to this as the partially-fixed initial condition.The final initial condition, of interest given that it greatly simplifies some of the analysis, is where all compartments are initiated from the joint stationary distribution for the system.We refer to this as the stationary initial condition and the system as a whole in this case as the stationary system.
3 Results and Discussion

Preliminaries
The multivariate Ornstein-Uhlenbeck process conditioned on the initial condition X 0 has exact solution [26,27] where and where ⊕ is the Kronecker sum.It follows directly that the stationary distribution, should it exist, is given by lim t→∞ We highlight that the non-stationary covariance matrix (eq.(5b)) does not depend on the initial condition X 0 and that the mean m(t) is an affine transformation of the initial condition X 0 .Therefore, for X 0 ∼ N (m 0 , Σ 0 ), we have that This expression reduces to the fixed initial condition for Σ 0 = 0, to the semi-fixed initial condition for Σ 0 = diag(σ 2 /(2θ), 0, . . ., 0), and to the stationary initial condition for Σ 0 = Σ ∞ .
The final result for the multivariate Ornstein-Uhlenbeck process that is relevant is the joint distribution of X t 1 ,t 2 ,... = X(t 1 ), X(t 2 ), . . .⊺ , which is multivariate normal with covariance matrix If t i ≫ 0 for all i such that Σ(t i ) = Σ ∞ , then eq. ( 8) corresponds to the joint stationary distribution of X t 1 ,t 2 ,... and, along with E(X(t i )) = µ, fully defines the system as a stationary Gaussian process.Furthermore, we can derive the distribution for all initial conditions (i.e., eqs.(5b), ( 6) and ( 7)) from the joint stationary distribution from marginalising or conditioning the joint stationary distribution accordingly (this is straightforward for the multivariate normal distribution, see [29]).

Quantifying smoothing in linear multicompartment processes
The structure of S, whereby noise enters the system only through the first compartment independently of other compartments, results in a simpler form of the stationary covariance matrix, given by eq. ( 6) and which we now denote simply by Σ ∞ , compared to the standard multivariate Ornstein-Uhlenbeck process.In particular, Σ ∞ depends only upon the first column of (Θ ⊕ Θ) −1 .In the supplementary material, we provide a full derivation for analytical expressions for Σ ∞ in two cases: the first where θ = k, and the second where both θ and k are allowed to vary freely.In this section, we summarise and discuss the main results.
For θ = k, elements of the symmetric matrix Σ ∞ are given by the recurrence relation subject to the boundary conditions The recurrence relation yields for covariances relating to the compartments themselves (i.e., i, j = 2, 3, . . .).Thus, the stationary variance of compartment ν ≥ 1 is given by where we have applied Stirling's approximation [30] to derive an asymptotic expression for the large compartment number, ν ≫ 1, limit.In fig.2a, we compare both the exact and asymptotic expressions for the stationary variance to a numerical approximation produced through repeated simulation of the SDE.Even for ν = 1, the asymptotic expression produces excellent agreement with simulation results.
Relaxing the restriction that θ = k yields a closed form solution for Σ ∞ , which simplifies along the diagonal to yield where B(•, •) and B(•, •, •) refer to the Beta function and incomplete Beta function, respectively.
We show both analytical and simulation results for σ 2 ν in this more general case in fig.2a.Taken all together, the results in fig.2a show that the variance dissipates as the compartment number increases.However, this happens relatively slowly: the analytical expression for θ = k provides a rate of decay of order ν −1/2 .Importantly, the expression in eq. ( 12) indicates that the variance does, indeed, tend to zero as ν → ∞.To the best of our knowledge it is not possible to derive a similar expression for general θ, however we provide in the supplementary material a simple proof that σ 2 ν+1 < σ 2 ν for θ > 0 to show that the variance is a strictly decreasing function and in fact tends to zero for large compartment numbers.From this, we also gain insight into the asymptotic behaviour of the solution for small or large θ.As θ → ∞, σ 2 0 → 0, and so the Shown are mean ± std from numerical simulations constructed from 10 replicates of 100 simulations (grey), and the corresponding analytical solution (colour).(b) Simulated and analytical ACFs for compartment ν = 3. (c) Second derivative of the ACF at ℓ = 0 for θ = k = 1 calculated directly from eq. ( 14) (diamonds), and using the analytical expression in eq. ( 16) (solid).(d-f) Analytical ACFs for I(t) (grey dashed) and X ν (t) (colour of increasing brightness for ν = 1, 2, . . ., 6).Unless otherwise stated, the other parameters are fixed at θ = µ = k = 1, and σ = 0.5.system becomes fully deterministic.In the problem itself, this represents a large mean reversion strength, so that effectively I(t) ≡ µ ∀ t.For θ = 0, the input function becomes purely Brownian motion and the stationary solution does not exist.
While compartment variances and covariances are statistics relating to the process at a single time point, the autocorrelation function (ACF) provides insight into how smooth the resultant time-series is.By deriving an analytical expression for e −Θ , we can obtain an analytical expression for the autocorrelation function (ACF).While this is relatively straightforward for the θ = k case, it is, however, significantly more involved in the more general case, for which the expression obtained is no more helpful that the analytical expression for the autocorrelation function involving the calculation of e −Θ directly (supplementary material).For θ = k we obtain where 1 F 1 ( • • • ) is the confluent hypergeometric function.Equation ( 14) can be expressed as where coefficients c i,ν depend on i and ν, which yields a simple expression for ν = 1.We also obtain the scaling of the autocorrelation as a function of ν by calculating the curvature of the ACF at ℓ = 0, where ′ indicates a derivative with respect to the lag, ℓ.In fig.2b, we compare analytical expressions for the ACF to those obtained through simulation, and in fig.2d-f we show the ACF for the system in fig. 1 for θ = 0.5, 1 and 2. In fig.2c we compare the analytical expression for the ACF curvature (eq.( 16)) to that calculated directly from eq. ( 14) using numerical methods.
As expected, we see qualitatively from the results in fig. 2 that further compartments remain correlated for longer.Considered alone, the stationary process X ν (t) is itself an, albeit non-Markovian, Gaussian process, fully defined by its variance and autocorrelation function.Therefore, should we normalise each compartment by its respective standard deviation, X ν (t)/σ ν (t), the properties of the resultant process are encoded entirely in the ACF.Given the system in eq. ( 1), we expect X ν (t) to be ν-times differentiable (the input, I(t) = X 0 (t) is nowhere differentiable) and therefore expect that further compartments will be smoother.For ν ≥ 1, we can see such smoothing directly from the ACF curvature in eq. ( 16).For a small increment, ℓ ≪ 1, where gives the dilation factor.Should we take ν 1 = 1, then X ν (t) varies a factor of ω 1,ν = √ 2ν − 1 more slowly than X 1 (t).

Intrinsic noise
By modelling the internal dynamics using a system of deterministic ODEs, we have implicitly assumed that the dimensional size of each compartment is sufficiently large relative to the input that we can neglect intrinsic noise.In this section, we relax this assumption and apply the linear noise approximation [31] (also known as the system-size expansion) to study the relative contribution to the stationary variance from both fluctuations in the input and from intrinsic noise arising for intermediate system sizes.We assume that Xν (t) = V X ν (t) corresponds to a dimensional concentration of material, where Xν (t) ∼ O(V ) such that V is the system size.
The governing equations for the internal compartments become The system can still be viewed as multivariate Ornstein-Uhlenbeck process, although S (eqs.(2) and ( 3)) is no longer a single-element matrix but has elements on both the main and lower diagonal.
As fluctuations driven by the Wiener processes W ν for ν ≥ 1 are independent of fluctuations in the input, we can decompose the stationary variance into contributions each from the input signal and intrinsic noise.The full derivation of the stationary variance for systems with finite Stationary Std . Adjustment to stationary variance due to intrinsic noise.Stationary standard deviation as a function of compartment number for (coloured curves) systems of various system sizes.The V → ∞ limit (grey solid) corresponds to eq. ( 13).Also shown are mean ± std from numerical simulations constructed from 20 replicates of 200 simulations (grey), and asymptotic approximations constructed using Stirling's formula (black dashed).Other parameters are fixed at θ = µ = k = 1, and σ = 0.5.
V are given as supplementary material.In summary, we stationary variance is now given by Eq. ( 13) Of particular note is that the contribution from intrinsic noise is both independent of k and bounded below by a minimum contribution of 1/(2V ) for ν = 1.In fig.3, we compare simulation results to both eq.( 19) and the large system-size case, V → ∞, which we denote as before by For sufficiently large ν, we can again apply Stirling's approximation [30] to see that the contribution from intrinsic noise behaves like Thus, while σ 2 ν vanishes, σ 2 intr,ν tends towards 1/V as ν → ∞.Moreover both these limits are approached at rates of order ν −1/2 .For all compartments, the contribution to the stationary variance from intrinsic noise is O(V −1 ), which compares (for the θ = k case) to the contribution from the input of O(σ 2 θ −3 ν −1/2 ).Therefore, the model that neglects intrinsic noise is valid not only for systems with sufficiently large system sizes, but also for systems with low compartment numbers where σ 2 θ −3 ≫ V −1 .

First passage time (FPT)
Motivated by the viral replication problem, we are now interested in studying the FPT distribution of individual compartments.That is, the time at which X ν (t) first crosses the threshold value X ν (t) = a from below, for a > X ν (0).To effectively compare FPT distributions between compartments, we scale the threshold by the associated stationary standard deviation of the relevant compartment such that a = µ + ãσ ν , where σ ν is given by eq. ( 12) and ã is specified.
Formally, we define the FPT by τ = inf{t : X ν (t) ≥ a} and its associated probability density (b,e) Distribution function for the FPT constructed from (colour) 1000 realisations of the SDE and (dashed black) a finite difference solution to eq. ( 21).(c,f) Mean, 2.5% quantile, and 97.5% quantile for the FPT distribution constructed from (grey) 1000 realisations of the SDE and (black) a finite difference solution to eq. ( 21).Shown in red-dashed is an approximation to the mean FPT constructed by scaling the FPT for ν = 1 based on matching the second-derivative of the autocorrelation function (eq.( 16)).The barrier for each compartment is located at a = ãσ ν where ã = 1.Other parameters are fixed at θ = µ = k = 1, and σ = 0.5.
and distribution functions by f (τ ) and F (τ ), respectively.We focus on fixed and partially-fixed initial conditions where we ensure that X ν (0) < a, demonstrated in fig.4a,d, respectively, for ν = 3 and ã = 1.
It is not generally possible to derive an analytical expression for f (t), nor to formulate a closed-form integral equation that can be solved numerically.The only general way to solve for f (t) is through a numerical solution to the Fokker-Planck equation, a ν + 1 dimensional partial differential equation.However, following the procedure for Markovian Gaussian processes [17,18] we are able to formulate a reasonable approximation for ν = 1, although we revert to simu lating the FPT for the more general case.
Denote by p ν (x, t) the density of the random variable X ν (t), and by p ν (x, t, τ ) the joint density with the first passage time.By marginalising p ν (x, t, τ ) with respect to τ we have that where K(t, τ ) = p ν (a, t|τ ) is the density of X(t) given the first passage time τ .The upper limit of t in the second integral arises since p ν (a, t|τ ) = 0 for all t < τ ; in other words, should a passage not have occurred by time t, then X ν (t) < a and so the particle cannot be at location 21) is a Volterra equation of the first kind and, while generally difficult, can be solved numerically provided K(t, τ ) can be computed.As I(t) = X 0 (t) is Markovian, p 0 (a, t|τ ) = p 0 (a, t|a, τ ) is readily available and for certain parameter combinations, eq. ( 21) can be solved analytically to give the FPT density of the Ornstein-Uhlenbeck process [32,33].
For ν > 0 the conditional probability p ν (a, t|τ ) cannot be calculated exactly; further we find that the standard approach of approximating p ν (a, t|τ ) ≈ p ν (a, t|a, τ ) [17] provides relatively poor results.To obtain a more accurate approximation, we note that the full state process X(t) is Markovian, and that X ν (s) = a and X ′ ν (s) > 0 if and only if s is a passage time.By eq. ( 1), X ′ ν (s) is a linear combination of other states, and so the random variable X ν (s), X ν (s) is Markovian with a multivariate normal distribution.Thus, we approximate Note that eq. ( 22) is not exact despite the full state being Markovian as we have not conditioned on a point, but rather the range X ′ (τ ) > 0; while the distribution of X ′ (s) is normal, the distribution of X ′ (τ ) is not necessarily so.In future, additional approximations based on socalled FPT functionals could potentially be constructed [34].We find that a numerical solution to eqs. ( 21) and ( 22) gives a reasonable approximation to f (τ ) for ν = 1 (fig.4b).
Results in fig.4b,e show the FPT distribution function, F (τ ), for both the fully-and partially-fixed initial conditions for θ = k = 1, respectively.The coloured curves are produced from 1000 realisations of the SDE model, and the black curves from a numerical solution to eq. ( 21).An interpretation of S(τ ) = 1 − F (τ ) is that of the survival probability: should a passage indicate system failure, S(τ ) represents the probability that a system is functional at time τ .For the virus-cell lysis problem, we interpret this as the probability that cell lysis has not occured, and viral production is ongoing.Visual inspection of the results in fig. 4 reveal little difference between both initial conditions, particularly for larger compartment numbers.
This observation is unsurprising upon comparison between the magnitude of the mean FPT, E(τ ) ∼ O (10), and the largest eigenvalue of −Θ, λ = −1, demonstrating that the influence of the initial condition decays like exp(−t) (eq.(5a)), much faster than the mean FPT.
The most obvious result from fig. 4b,e, as one might expect from the analysis of compartment smoothing in the previous section, is that the FPT is generally larger for further compartments; accounting for differences in the stationary variance by comparing compartments across the same value of ã indicates that further compartments can be thought to be more robust to external noise.Not only does the expected FPT increase (equal to the area under the survival function S(τ )), but so too do the lower quantiles, evidenced by the time taken for the distribution function F (τ ) to visually become non-zero.We investigate these qualitative observations further in fig.4c,f, by calculating the mean, 2.5% and 97.5% quantiles for the FPT for each compartment.
Aside from an overall increase in the FPT, there is a notable increase in the inter-quantile range or variance for larger compartment numbers.
By normalising the location of the threshold by the compartment standard deviation, the FPT is almost entirely a function of the ACF for each compartment.The compartment variance still plays a role for the initial conditions considered in fig.4, as the relative distance between the initial condition X ν (0) = 0 and the threshold depends on σ ν .The settling phase, however, occurs relatively quickly: for k = 1, the system settles to equilibrium like exp(−t), much faster than the typical FPT.Thus, the compartment smoothing effect characterised by ρ ν (ℓ), which provides the temporal scale, is the primary factor that drives increases in FPT for further compartments.In fig.4c,f, we show that the small-lag ACF dilation factor, given analytically by eq. ( 17), provides an excellent match to numerical results for the mean first passage time, confirming qualitatively that, similarly to the compartment variance, the FPT scales like √ ν.
Additional results (supplementary material), constructed for various values of k and hence various ACFs, demonstrate strong correlation (Spearman correlation of −0.987) between the ACF curvature and the expected FPT.

Production before system failure
While for many systems the FPT, τ , may itself be of primary interest, for others it may be a so-called FPT functional that is important [34].For the virus-cell lysis problem, for example, of primary interest may be the total amount of virus produced by the system prior to cell lysis: a viral genotype that maximises per-host-cell virion production may have a fitness advantage over others.Also seen in the virus literature is that this so-called "burst-size" declines with cell size [35], potentially linking the structure of a particular viral replication network to an optimal cell size and, by extension, an optimal cell type.
We define cumulative production as the total amount of material to pass out of the system, A(t) = t 0 kX ν (t) dt, and investigate A(τ ) through simulation in fig. 5 for ν = 6 and θ = k = 1.Results in fig.5a show that A(τ ) is highly correlated with τ .We expect this, as both initial settling and ACF decay occur much faster on average than τ .Thus, for the linear system considered in this work, we hypothesise that a genotype that maximises the expected FPT could be considered equivalent to one that maximises the production prior to lysis.
In fig.5b,c we perform a parameter sweep to determine the distribution of A(τ ) as a function of the threshold location, ã, and the compartment transfer rate, k.Clearly, for thresholds that are larger, and crossed more infrequently, we see an increase in A(τ ).We view the threshold location as a feature of the host-system: for the virus-cell lysis problem, this is a biological feature of the host cell, and not a feature of the viral genome.Of direct interest is the relationship between k and A(τ ).The results in fig.5c show that systems with large compartment transfer rates have lower cumulative production than those that operate slowly.While we cannot interpret the expression for the ACF curvature (eq.( 16)) exactly for θ ̸ = k, the expression does suggest that the ACF curvature is proportional to k 2 , thus increasing the rate at which material travels through the system is detrimental to robustness.
In the virus-cell lysis problem, these results appear to suggest that production can be maximised should the system tune itself to operate slowly (i.e., small values of k).However, this observation ignores potential tradeoffs induced by immune-response mechanisms in the host [36]: while production in an isolated or experimental system might be maximised by slowly material progression, operating quickly carries advantages of lysis occuring due to material production, and not through an immune response of the host organism.Analysis of more complicated, non-linear, compartment processes may also yield non-monotonic relationships between fitness and speed other parameters; such analysis is, however, beyond the scope of the present work.

Non-local feedback to alter system robustness
Another way in which systems can potentially increase their robustness to input noise is through non-local feedback or feedforward loops.Such a phenomenon, where progression through the virus life cycle is not unidirectional, is with evidence in the virus literature: for example, the complex network-like replication cycle seem in the human ademovirus [37].Hepatitus B, an enveloped DNS virus, recycles new viral DNA that is awaiting repackaging back into the nucleus [38,39]: analogous, in our simple linear model, to a feedback from the final to an early compartment.A similar late to early-stage feedback is seen in positive-strand RNA viruses in which newly replicated RNA strands are either encapsidated or reutilised in translation and replication [40].More generally, positive-and negative-strand RNA viruses have very similar replication structures and an identical feedback mechanism, however are in part distinguished by their feedforward structures [41,42].Less concretely, replication of the measles virus comprises several feedback and feedforward mechanisms, particularly during the transition and polymerase stages [43].
In this section, we investigate the relative change to the stationary variance and ACF curvature in the final compartment of a system with an additional transfer from compartment n to compartment m of magnitude ε (fig.6a).A transfer with m < n is considered a feedback, with m > n a feedforward.Since the within system dynamics are deterministic, a transfer with m = n has no net effect on the dynamics.
The results in fig.6b,c show that the output variance can be reduced, potentially significantly, through a feedback.For the system considered, a feedback from the final to the first compartment has the largest effect, yielding a reduction of over 30% to the stationary standard deviation.In fig.6d,e we show a similar affect on the curvature of the ACF.Interestingly, feedbacks from the second last compartment to the first compartment, or the last to the second, yield the greatest reduction in the curvature of the ACF.We also observe a non-linear relationship between the magnitude of the feedback and the ACF curvature.For example, a -5% -2.5% 0 +2.5% +5% -20% 0 +20% X 1 X m X n ...  small (ε = 0.1) feedforward rate from the first to the last compartment yields a reduction in the ACF curvature whereas a large (ε = 1) feedforward rate yields an increase in the curvature.

Continuum limit
While the focus of the present work is on multicompartment processes, a natural extension is to investigate smoothing in processes that occur on a continuum (for example, where compartment number, ν, is a continuous measure of how far a particle has progressed through a system).
Therefore, we consider a refinement of the discrete process by dividing each compartment into 1/∆ subcompartments, each of width ∆ (fig.7a).To maintain the effective total time a particle spends in the system, we assume that the subcompartment transfer rate becomes k = k/∆ such that X 0 (t) = I(t), where ν = (i−1)∆.Note that to formulate the input as a boundary condition, we have modified the original system such that the input transfers to the first compartment at rate k (i.e., X 1 (t) now experiences an input of kI(t) compared to the input of I(t) in eq. ( 1)).This formulation is equivalent to the original formulation for k = 1.
We denote x(ν, t) = X i (t), and take ∆ → 0 to yield an advection equation with Dirichlet X 1 X 2 X n ... boundary condition with exact solution As an advection equation, this continuum analogue of the multicompartment system corresponds to exact (undamped) transport of material through the system.We conclude, therefore, that the smoothing we see is a uniquely discrete effect.This conclusion becomes obvious should each compartment be viewed as a well-mixed segment of "length" ∆ in ν-space.Transfer between each compartment represents flow across the left boundary.As the "length" of each compartment becomes smaller, the left boundary approaches the right, and thus the finite difference between successive compartments diminishes.Thus, smoothing can be viewed as a consequence of within-compartment mixing, a feature of the discrete system that vanishes as the size of each compartment tends to zero.
We perform a numerical experiment in fig.7b-d and simulate three systems subject to an identical input, with identical total length n = 10, and with compartment spacing reducing from ∆ = 1 to ∆ = 0.01.In effect, the solution to the discrete system corresponds exactly to a forward difference approximation to that of the advection equation (eq.( 24)).Smoothing in both the autocorrelation and variance is evident in all cases with finite ∆, however, it diminishes significantly for ∆ = 0.01.
The advection equation given in eq. ( 24) gives little insight into the behaviour as ∆ → 0.
To derive a continuum limit approximation that captures the smoothing effects in systems with small but finite ∆, we apply the method of multiple scales and choose a 'slow' scale of This scaling yields or, in the original variables, Full details relating to the derivation of this high-order continuum limit approximation are given as supplementary material.
In fig.7e, we provide a numerical solution to eq. ( 27) for ∆ = 1, showing that this second continuum approximation captures the smoothing effect seen in the discrete model.The inclusion of the diffusion term with coefficient O(∆) demonstrates that, while the problem is singular, smoothing is always present in the discrete system.This observation is, in fact, a well known feature of the discrete system when viewed as an approximate numerical solution to the advection equation with an upwind spatial scheme.Such a scheme can be simply seen to introduce artificial diffusion-effectively, smoothing-with form identical to that given by the multiple scales approach in eq. ( 27).

Conclusion
Multicompartment processes are ubiquitous in biology; from linear progression through the cell cycle, to phage replication in bacteria and the propagation viruses by hijacked cellular machinery.Our analysis demonstrates that even a fundamental linear multicompartment structure provides potential advantages and benefits to the systems that employ them.These results parallel filters and autoregressive models in control theory that allow engineers to control and exploit systems subject to noisy input [25,44,45].
Most notable is the effect of such systems to smooth and provide an additional degree of control over external noise, consequentially increasing resilience and robustness.The inclusion of feedback and feedforward loops can enhance this effect, providing systems with additional degrees of control and contributing to so-called perfect adaptation [45,46].Such loops provide a potential explanation for out-of-order progression in some systems, for example, whereby viral replication does not occur as a unidirectional process [2].Our work demonstrates that feedback loops could yield a fitness advantage through more favourable statistical properties of viral load compared with perfect progression through the replication cycle.Such results potentially explain the complexity in the replication network structure seen in some viruses; for example, in the human adenovirus [37].Indirectly, these loops (and by extension, more complicated network structures) provide systems with the ability to tune the first passage time distribution, potentially yielding an optimal lysis time [13,47].While we restrict our analysis to a single nonlocal feedback or feedforward, future work may study more general optimal network structures informed by more specific biological problems.
Analysis of a linear model system, subject to Ornstein-Uhlenbeck-type noise, allows us to present closed-form expressions for key statistics, providing a fundamental understanding that would be otherwise unavailable for more complicated systems.We reveal that noise dissipates in systems with unidirectional flow, eventually vanishing in an infinite-compartment system.However, we show that this effect is tied to a finite flow rate: scaling to yield continuous flow through a continuum limit approximation reveals that smoothing is a discrete effect, caused by within compartment mixing.While many simple results are only available for the constant flow rate assumption, arbitrarily connected linear systems yield a multidimensional Ornstein-Uhlenbeck process with statistical properties computable semi-analytically through the exponentiation of a connectivity matrix.
Importantly, we lay the foundation for future work to explore the properties of more general multi-compartment systems subject to external noise.For instance, study of the interaction between the external timescales (i.e., autocorrelation of the external input) and the internal timescales (i.e., progression through the system) is highly relevant to specific biological problems: in the virus-cell lysis problem, this interaction is relevant for immune system detection and, therefore, infection clearance.Aside from external noise, other stochastic features, including both between-cell and between-virion heterogeneity and fluctuations in the replication process itself are known to play an important role in within-host virus replication [3,10,48,49].Despite these observations, the study of multicompartment problems with the stochastic mathematical models requisite to capture important features is presently scarce, albeit a rich area for mathematical and biological insight.

Figure 1 .
Figure 1.Multicompartment model of viral replication subject to stochastic input.(a) The general viral replication cycle.Virions enter a cell and progress through a multistage process, before accumulating following repackaging.(b)We study a linear analogue of the virus problem, namely a system subject to random input, I(t), modelled as an Ornstein-Uhlenbeck process with mean µ, reversion strength θ, and noise magnitude σ.A realisation of the input is shown in (c).Material then progresses through n compartments at constant rate k.(d) A realisation of the system in (b), where X ν (t) models the concentration of material in each compartment ν = 1, 2, . . ., 6.It is evident that passage through the compartments has a smoothing effect.

Figure 4 .
Figure 4. First passage time distributions for linear multicompartment model.(a,d) Realisations of a three compartment system initiated using (a) the fixed initial condition and (d) the partially-fixed initial condition.Solutions are terminated at t = τ : X 3 (τ ) > a, yielding τ as the FPT.(b,e)Distribution function for the FPT constructed from (colour) 1000 realisations of the SDE and (dashed black) a finite difference solution to eq. (21).(c,f) Mean, 2.5% quantile, and 97.5% quantile for the FPT distribution constructed from (grey) 1000 realisations of the SDE and (black) a finite difference solution to eq. (21).Shown in red-dashed is an approximation to the mean FPT constructed by scaling the FPT for ν = 1 based on matching the second-derivative of the autocorrelation function (eq.(16)).The barrier for each compartment is located at a = ãσ ν where ã = 1.Other parameters are fixed at θ = µ = k = 1, and σ = 0.5.

Figure 5 .
Figure 5. Cumulative production at the first passage time.We investigate the cumulative production at the first passage time, denoted by A(τ ) = τ 0 kX 6 (t) dt, for a six compartment system.(a) Relationship between τ and A(τ ) based on 1000 replicates of the SDE with ã = k = 1.We show both a scatter plot of the joint density and kernel density estimates constructed from the marginals.(b,c) Mean, 2.5% quantile, and 97.5% quantile of A(τ ) constructed from 1000 replicates of the SDE as (b) ã is varied, and (c) k is varied.

Figure 6 .
Figure 6.System with non-local feedback.We investigate the effect of a non-local feedback (or feedforward) of magnitude ε from compartment n to compartment m on the stationary standard deviation of X 6 (t), denoted σ 6 , and the magnitude of the ACF curvature, denoted |ρ ′′ 6 (0)|.In all cases, a reduction in each statistic corresponds to a potentially more robust system.

Figure 7 .
Figure 7. Material transport through a system approaching the continuum limit.(a) Compartments in the discrete system are divided into subcompartments of "length" ∆, while the "transfer density" is kept fixed.(b-d) Solutions of the discrete system for decreasing ∆.All systems are subject to identical input I(t) (grey).In (d), the shifted input I(t − 10) is shown at ν = 10 for comparison with the transported concentrations.Other parameters are fixed θ = µ = k = 1, and σ = 0.5.(e) Solution to the continuum limit approximation (eq.(27)) for ∆ = 1.