Renewal reward perspective on linear switching diffusion systems
Maria-Veronica Ciocanel, John Fricks, Peter R. Kramer, Scott A. McKinley
NNoname manuscript No. (will be inserted by the editor)
Renewal reward perspective on linear switchingdiffusion systems
Maria-Veronica Ciocanel · John Fricks · Peter R. Kramer · Scott A. McKinley
Received: date / Accepted: date
Abstract
In many biological systems, the movement of individual agents ofinterest is commonly characterized as having multiple qualitatively distinctbehaviors that arise from various biophysical states. This is true for vesiclesin intracellular transport, micro-organisms like bacteria, or animals movingwithin and responding to their environment. For example, in cells the move-ment of vesicles, organelles and other intracellular cargo are affected by theirbinding to and unbinding from cytoskeletal filaments such as microtubulesthrough molecular motor proteins. A typical goal of theoretical or numericalanalysis of models of such systems is to investigate the effective transport prop-erties and their dependence on model parameters. While the effective velocityof particles undergoing switching diffusion dynamics is often easily character-ized in terms of the long-time fraction of time that particles spend in eachstate, the calculation of the effective diffusivity is more complicated becauseit cannot be expressed simply in terms of a statistical average of the parti-cle transport state at one moment of time. However, it is common that thesesystems are regenerative, in the sense that they can be decomposed into inde-pendent cycles marked by returns to a base state. Using decompositions of this
MVC is supported by The Ohio State University President’s Postdoctoral Scholars Programand by the Mathematical Biosciences Institute at The Ohio State University through NSFDMS-1440386. JF, PRK, and SAM are supported by NIH R01GM122082-01.Maria-Veronica CiocanelMathematical Biosciences Institute, The Ohio State UniversityTel.: +614-688-3334E-mail: [email protected] FricksSchool of Mathematical and Statistical Sciences, Arizona State UniversityPeter R. KramerDepartment of Mathematical Sciences, Rensselaer Polytechnic InstituteScott A. McKinleyDepartment of Mathematics, Tulane University a r X i v : . [ q - b i o . S C ] N ov Maria-Veronica Ciocanel et al. kind, we calculate effective transport properties by computing the momentsof the dynamics within each cycle and then applying renewal-reward theory.This method provides a useful alternative large-time analysis to direct homog-enization for linear advection-reaction-diffusion partial differential equationmodels. Moreover, it applies to a general class of semi-Markov processes andcertain stochastic differential equations that arise in models of intracellulartransport. Applications of the proposed renewal reward framework are illus-trated for several case studies such as mRNA transport in developing oocytesand processive cargo movement by teams of molecular motor proteins.
Keywords renewal rewards · intracellular transport · processive motortransport Mathematics Subject Classification (2000) · · · Microscale biological agents frequently change biophysical state, which resultsin significant changes in their movement behavior. Intracellular cargo, for ex-ample, switches among active transport, diffusive transport, and paused states,each resulting from different mechanochemical configurations of the cargo,cytoskeletal filaments, and the molecular motors that bind them (Hancock,2014; Bressloff and Newby, 2013). Models for this behavior can be either de-terministic (typically partial differential equations, PDEs) or stochastic (oftencontinuous-time Markov chains, CTMCs, or stochastic differential equations,SDEs) depending on whether the investigation focuses on population proper-ties (deterministic methods) or individual paths (stochastic methods). Eachstate is commonly characterized in terms of a mean velocity, fluctuations aboutthe mean velocity, and a distribution of time spent in the state, sometimes butnot always determined by classical reaction rate theory. Explicit solutions forthese models are rarely available, so asymptotic or numerical methods areoften deployed to investigate and characterize the model’s predictions. Thestudy of deterministic models often relies on numerical simulation using PDEintegration methods (Wang et al, 2003; Cox and Matthews, 2002; Trong et al,2015), while stochastic models are simulated with Monte Carlo/Gillespie al-gorithms (M¨uller et al, 2008; Kunwar and Mogilner, 2010; M¨uller et al, 2010;Allard et al, 2019) to generate individual trajectories that are then analyzedstatistically. However, these computations can be quite costly, especially whenone wants to understand how bulk transport properties (like effective velocityor diffusivity) depend on individual model parameters. When possible, asymp-totic analysis allows for explicit approximation of transport properties, whichcan validate, complement, or even replace numerical simulations (Reed et al,1990; Brooks, 1999; Pavliotis, 2005; Pavliotis and Stuart, 2008; Popovic et al,2011; McKinley et al, 2012; Bressloff and Xu, 2015; Ciocanel et al, 2017).The long-term effective velocity of state-switching particles is often straight-forward to compute, usually obtained by calculating the fraction of time spent enewal reward perspective on linear switching diffusion systems 3 in each state and correspondingly averaging the associated state velocities. Onthe other hand, this weighted average technique is not valid when calculatinga particle’s effective diffusivity, since this quantity cannot be simply expressedin terms of the statistics of the particle transport at a single moment of time.Rather the effective diffusivity depends on temporal correlations that existin displacements, including across changes in biophysical state. Generalizingsome previous work (Brooks, 1999; Hughes et al, 2011; Krishnan and Epure-anu, 2011; Hughes et al, 2012; Ciocanel et al, 2017), we consider this problemof computing effective diffusivity for a class of state-switching particle modelsthat can be expressed in a framework where the sequence of states are givenby a Markov chain, but the times spent in these states are not necessarilyexponentially distributed as in a continuous-time Markov chain. Since we as-sume that the state process Markov chain is positive recurrent, the particleposition process can be described as a regenerative increment process in asense defined by Serfozo (2009), for example. That is to say, we consider pro-cesses that almost surely return to some base state at a sequence of (random)regeneration times such that the dynamics after a regeneration time are in-dependent from those that occur before. As a result, we can decompose theprocess into what we refer to as cycles , in which the particle starts in a basestate, undergoes one or more state changes, and then revisits the base stateagain. The dynamics within each cycle are independent of other cycles andwe can use the renewal-reward theorem to perform asymptotic calculationsby viewing the total displacement within each cycle as its reward and view-ing the cycle durations as times between regenerations. An early applicationof the idea of computing effective particle velocity and diffusivity by decom-position and analysis of the dynamics in terms of independent cycles was tolarge enhancement of (non-switching) particle diffusion in a tilted periodicpotential (Reimann et al, 2002, 2001).Our primary motivating examples are related to intracellular transport.Some prominent recent investigations include the study of mRNA localizationin oocyte development (Zimyanin et al, 2008; Trong et al, 2015; Ciocanelet al, 2017), cell polarization in the budding yeast (Bressloff and Xu, 2015),neurofilament transport along axons (Jung and Brown, 2009; Li et al, 2014),interactions of teams of molecular motor proteins (Klumpp and Lipowsky,2005; M¨uller et al, 2008; Kunwar and Mogilner, 2010; M¨uller et al, 2010),and sliding of parallel microtubules by teams of cooperative identical motors(Allard et al, 2019). Microtubule-based transport of cargo is typically mediatedby kinesin motors moving material to the cell periphery and by dynein motorscarrying it to the nucleus. Understanding population-scale behaviors, such asprotein localization, that arise from local motor interactions remains an openquestion. While multiple motor interactions are usually thought to be resolvedthrough a tug-of-war framework (M¨uller et al, 2008), it has been observed thatimportant predictions made by the tug-of-war framework are not consistentwith in vivo experimental observations (Kunwar et al, 2011; Hancock, 2014).The work presented in this paper can aid theoretical efforts to relate localmotor-cargo dynamics to predictions for large scale transport.
Maria-Veronica Ciocanel et al. ∂ u ( y, t ) ∂t = A T u − V ∂ y u + D∆ u . (1)We will actually think of u as an ( N + 1)-dimensional column vector (in-dexed from 0 to N ) of the concentrations of particle populations in differentdynamical states, which also obey the forward Kolmogorov equations witha different normalization. The dynamics are governed by matrices A, V, D ∈ R ( N +1) × ( N +1) , where V and D are diagonal matrices, with real constant di-agonal entries v , v , . . . , v N for V corresponding to the particle velocities ineach state, and positive real constant diagonal entries d , d , . . . , d N for D corresponding to the diffusion coefficients in each state. The matrix A is thetransition rate matrix of the associated finite state continuous-time recurrentMarkov chain (CTMC), J ( t ), which tracks the state of the particle at a giventime. That is to say, each off-diagonal entry a ij can be interpreted as the rateat which a particle in state i switches to state j . The diagonal entries of A arenon-positive and correspond to the total rate out of a given state. The rowsof A sum to zero. Assuming that the CTMC is irreducible, it follows that A admits a zero eigenvalue with algebraic and geometric multiplicity one, andthe corresponding normalized zero-eigenvector π is the stationary distributionof J ( t ).Either quasi-steady-state reduction (Bressloff and Newby, 2013) or ho-mogenization theory (Pavliotis and Stuart, 2008) can be used to reduce thecomplexity of the advection-diffusion-reaction system 1 to a scalar advection-diffusion equation of the form: ∂c ( y, t ) ∂t = v eff ∂ y c ( y, t ) + D eff ∆c ( y, t ) . (2)with constant effective velocity v eff and constant effective diffusivity D eff forthe particle concentration without regard to state c ( y, t ) = (cid:80) Ni =0 u i ( y, t ).Quasi-steady-state reduction assumes the stochastic switching dynamics oc-curs on a fast scale relative to the advection-diffusion dynamics, while homog-enization theory applies at sufficiently large space and time scales relative tothose characterizing the dynamical scales. These different asymptotic condi-tions give in general distinct results when the transport and switching rateshave explicit dependence on space, but when, as in the present case, they arespatially independent, the formulas for the effective transport coefficients co-incide. (This is because the time scale of advection/diffusion is linked purely tothe spatial scale, so the large spatial scale assumption of homogenization willperforce induce a time scale separation between the switching and transport enewal reward perspective on linear switching diffusion systems 5 dynamics, as assumed in quasi-steady-state reduction.) The effective velocityis computed by averaging the velocity in each state, weighted by the stationarydistribution of the particle states: v eff = v · π , (3)where v = ( v , v , . . . , v N ) T . The effective diffusivity is given, from an equiv-alent long-time effective dynamical description for intracellular transport de-rived by Lyapunov-Schmidt reduction on the low wavenumber asymptotics ofthe Fourier transform of Eq. 1 in Ciocanel et al (2017, 2018), by D eff = d · π − v · ( A T ) − ( v ◦ π − v eff π ) , (4)with v ◦ π = ( v (0) π (0) , v (1) π (1) , . . . , v ( N ) π ( N )) T (5)denoting the Hadamard product (componentwise multiplication of vectors).Here d = ( d , d , . . . , d N ) T , and A T is the restriction of A T to its rangeRan( A T ) (vectors orthogonal to (1 , , . . . , T ). Note that the operation in-volving the inverse of ( A T ) − is well-defined since A T is a full-rank matrixmapping Ran( A T ) to Ran( A T ), and its inversion in Eq. 4 is applied to a vectorin Ran( A T ). We remark that the homogenization formula is often written (Cio-ranescu and Donato, 1999; Pavliotis and Stuart, 2008) in an equivalent adjointform to Eq. 4, with a centering of the leading vector v → v − v eff (1 , , . . . , T that renders the formula indifferent to the choice of how to invert A T . Theterm d · π above reflects the contributions to the asymptotic diffusivity frompure diffusion, while the second term captures the interactions between theadvection and reaction terms.Applications of quasi-steady-state reduction to biophysical systems withstate switching and diffusion can be found in Newby and Bressloff (2010b,a);Bressloff and Newby (2011, 2013); Bressloff and Xu (2015). Homogenizationof Brownian motor models was conducted in Pavliotis (2005); Kramer et al(2010).1.2 Summary of method based on regeneration cyclesThese foregoing methods (Pavliotis and Stuart, 2008; Ciocanel et al, 2017;Bressloff and Newby, 2013) rely on the fully Markovian structure of the dy-namics, with the state-switching process in particular taking the form of acontinuous-time Markov chain with exponentially distributed state durations.In this work, we consider a generalized framework in which we require only thatthe sequence of states visited form a discrete-time recurrent Markov chain, butdo not require exponentially distributed state durations, so the state-switchingprocess J ( t ) need not be a continuous-time Markov chain. Moreover, we allowfor more general random spatial dynamics within a state that also need notbe fully Markovian. Our framework and method of analysis rather requiresonly a regenerative structure of the dynamics, with repeated returns to a base Maria-Veronica Ciocanel et al. state, at which moments the future randomness is decoupled from the previoushistory.We use renewal-reward theory and a functional central limit theorem toderive effective drift and diffusion for these more general switching diffusionsystems in terms of the analysis of a single regeneration cycle. The calculationframework also results in an expression for the expected run length of cargoundergoing switching diffusion. Our approach builds on previous applicationsof renewal-reward processes modeling motor-stepping and chemical cycles ofbead-motor assays (Krishnan and Epureanu, 2011; Hughes et al, 2011, 2012;Miles et al, 2018; Shtylla and Keener, 2015) and extends the technique toaccommodate more complex models with dynamics depending on the amountof time spent in the current state, as described in §
2. Given the renewal-rewardframework, the analysis of the model reduces to computing the correlatedspatial displacement and time duration of each cycle, which we study in § §
4, we show that our method of deriving effec-tive velocity and diffusivity agrees with predictions in Ciocanel et al (2017)arising from a Liapunov-Schmidt reduction approach equivalent to homoge-nization for partial differential equations describing mRNA concentrations asin (1). In §
5, we show that our method also agrees with previous theoreticaland numerical analyses of transport properties for cargo pulled by teams ofmolecular motors. In the case of tug-of-war dynamics, with cargo transportedby teams of opposite-directed motors, our framework provides predictions onthe dependence of effective diffusivity on the ratio of stall to detachment forceof the pulling motors. We also apply this method to a model accounting forincreased reattachment kinetics when motors are already attached to the cargoand show that teams of opposite-directed motors have lower effective velocitiesbut larger run lengths than teams consisting of the dominant motor only. Fi-nally, we show that our effective diffusivity calculation agrees with stochasticsimulations of sliding microtubule behavior driven by teams of bidirectionalmotors for a large range of load sensitivity. As the experimental data on mo-tor interactions develops rapidly, the framework proposed may prove usefulin analyzing novel models and in understanding the dependence of effectivetransport properties on model parameters.
The type of path we have in mind in this work is displayed in Figure 1, acontinuous, stochastic process that switches between several stereotypical be-haviors. Let the real-valued process { X ( t ) : t ≥ } be the time-dependentposition of a particle and let { J ( t ) : t ≥ } denote the time-dependent un-derlying (e.g., biophysical) state, taking values from the finite state space S = { , , , . . . , N } . Switches between the states take place at the randomtimes { t k : k ∈ N } and we use { J k : k ∈ N } to denote the state during the k th time interval [ t k − , t k ). We set t = 0 and J = J (0). We assume that the enewal reward perspective on linear switching diffusion systems 7 sequence of states { J k : k ∈ N } forms a time-homogeneous recurrent Markovchain with zero probability to remain in the same state on successive epochs.Given the state J k , the associated state duration t k − t k − and spatial dis-placement X ( t k ) − X ( t k − ) are conditionally independent of all other randomvariables in the model (but not necessarily of each other). Moreover, the con-ditional joint distribution of t k − t k − and X ( t k ) − X ( t k − ) given J k dependsonly on the value of J k and not separately on the index k . In other words, thedynamics of ( J ( t ) , X ( t )) have a statistically time-homogeneous character. Time t P o s i t i on X ( t ) t t T ⌧ t T Fig. 1
An example of the type of particle trajectory considered in this work. Intracellularcargo in the form of a vesicle experiences periods of active and diffusive transport. Thedashed vertical lines indicate random times { t k : k ∈ N } when there are switches in thebiophysical state. The base “renewal” state is free diffusion and the red dashed vertical linesmark moments { T k : k ∈ N } when the system enters the base state. We denote the timesspent at each time step by τ k , as detailed in § One general subclass of the processes considered can be expressed as fol-lows: The random times { t k } ∞ k =0 are generated by sampling t k − t k − indepen-dently from their conditional marginal distributions given the Markov chainstates J k , and then conditioned upon these random variables, the spatial pro-cess X ( t ) is governed by a stochastic differential equation with coefficientsdepending on the current state, the value of X upon entry into the currentstate, and the time since entry into the current state. That is, we express theconditional dynamics of X ( t ) as:X. ( t ) = ∞ (cid:88) k =1 [ t k − ,t k ) ( t ) (cid:16) α J k (cid:16) X ( t ) , X ( t k − ) , t − t k − (cid:17) t. + (cid:112) d J k W. ( t ) (cid:17) , (6)where α j : R × R + → R is a function that describes the drift (deterministiccomponent of the dynamics) of a particle while in the state j , and d j is thediffusivity in that state. (In general, the diffusive coefficients might also dependon the position of the particle and the recent history of the process, but we Maria-Veronica Ciocanel et al. restrict ourselves to memoryless, additive noise for this discussion.) In thefollowing paragraphs, we describe a few examples.First, we can consider the stochastic process associated with the PDE(1), for which we set the drift terms in (6) to be α j = v j , where v j is theconstant j th diagonal entry of the velocity matrix V in (1). The diffusioncoefficients would correspond to the entries of the diagonal matrix D in (1).The process would switch between states as a CTMC with rate matrix A ,which, through standard theory (Lawler, 1995, Sec. 3.2), induces a transitionprobability matrix for the sequence of states { J k : k ∈ N } : P ij = (cid:40) A ij ¯ λ i for i (cid:54) = j, i = j. where ¯ λ i ≡ (cid:80) j ∈ S \{ i } is the total transition rate out of state i . The stateduration in state i would be exponentially distributed with mean ¯ λ − i .For a second example, we consider the process depicted in Figure 1. Thisprocess could describe a model for intracellular transport, with transport statesdriven by processive molecular motors. There are four states: a forward pro-cessing state, a backward processing state, and a stalled state – each of whichare characterized by having a drift with speed v j plus Ornstein-Uhlenbeck typefluctuations (as described for example in Smith and McKinley (2018)) – anda freely diffusing state where the drift term equals zero. That is, α = 0 forthe freely diffusing state and for j > α j ( y, y , t ) = − κγ (cid:0) y − ( v j t + y ) (cid:1) , (7)where κ is a spring constant, γ is the viscous drag and v j is the velocityassociated with the j th state. The term ( v j t + y ) indicates the theoreticalposition of a processive molecular motor that is simultaneously bound to theparticle and to a microtubule. Remark 1
We note that there are at least two ways that the process X ( t )can be considered to be non-Markovian and still fall within the set of modelsto which our results apply. The first, which is captured by the drift term(7), is that the process X ( t ) has memory in the sense that resolving X ( t )on any interval in ( t k , t k +1 ) depends on the value X ( t k ). A second allowablenon-Markovian dynamic can be obtained by choosing the state duration times t k − t k − given state J k to have a non-exponential distribution. As long as thestochastic process of states { J k } is a time-homogeneous, positive recurrentMarkov chain, the technique we present will apply.In § enewal reward perspective on linear switching diffusion systems 9 by Brooks (1999) in which each state is characterized by a fixed velocity pa-rameter α j = v j with the state diffusivity d j set to zero. In some of the otherexamples, fluctuations about the mean are included and would contribute tothe long-term diffusivity as a result. Of course, fluctuations are always presentin these dynamics (sometimes due to variability in the motor stepping, some-times due to fluctuations in cargo position). There are natural ways to addthese considerations to the models in § regeneration times { T n } . In what follows, we will view the con-secutive spatial displacements and time durations of the regenerative cycles tobe the rewards and cycle durations of a classical renewal-reward process (Cox,1962). Because the cycle statistics are independent and identically distributedafter the first regeneration time T , we define (in the sense of distribution)random variables for a generic cycle n ≥ ∆X D = X ( T n ) − X ( T n − ); ∆T D = T n − T n − ; and M D = sup t ∈ [ T n − ,T n ] | X ( t ) − X ( T n − ) | . (8)We rely on the functional central limit theorem (FCLT) presented in Ser-fozo (2009) for our asymptotic results. To this end, we define the quantities µ := E ( ∆T ); a := E ( ∆X ) E ( ∆T ) ; and σ := Var( ∆X − a∆T ) . (9)As in previous work on molecular motor systems (Hughes et al, 2011, 2012),the FCLT justifies defining the effective velocity and effective diffusivity of theprocess X ( t ) in terms of properties of the regenerative increments as follows: v eff := lim t →∞ t X ( t ) = a = E ( ∆X ) E ( ∆T ) ; (10) D eff := lim t →∞ t Var( X ( t )) (11)= σ µ = 12 E ( ∆T ) (cid:16) Var( ∆X ) + v Var( ∆T ) − v eff Cov( ∆X, ∆T ) (cid:17) . In more technically precise terms, the FCLT states: For r ∈ Z + , define Y r ( t ) :=( X ( rt ) − art ) / ( σ (cid:112) r/µ ) . If a, µ, σ, E ( M ), and E (( ∆T ) ) are all finite, thenlim r →∞ Y r = B in distribution for t ∈ [0 , { B : t ∈ [0 , } is astandard Brownian motion (Whitt, 2002).2.2 Notation for events within each regeneration cycleThe mathematical analysis in § § n th cycle by η ( n ) := min { k ≥ J k + K n − +1 = 0 } , where K = 0 and K n = (cid:80) ni =1 η ( i ) . We will let τ ( n ) k = t K n − + k − t K n − + k − denote the times spent in each step of the n th cycle, and ξ ( n ) k = X ( t K n − + k ) − X ( t K n − + k − ) denote the corresponding spatial displacements. The total time ∆T and displacement ∆X accrued in a cycle n ≥ ∆T := η ( n ) (cid:88) k =1 τ ( n ) k and ∆X := η ( n ) (cid:88) k =1 ξ ( n ) k . (12)In what follows, we drop the superscript denoting the index n of the cycle, sincethe cycles have statistically independent and identically distributed behaviorfor n ≥
2. We will decompose each cycle into what is accrued during the firststep ( τ and ξ ) associated with the visit to the base state, and what accruesin all subsequent steps in the cycle, which we label ∆ ˜ T := ∆T − τ and ∆ ˜ X := ∆X − ξ . (13)For each state j ∈ S of the underlying Markov chain, let { τ k ( j ) , ξ k ( j ) } ∞ k =1 be a sequence of iid pairs of random variables drawn from the conditionaljoint distribution of durations and displacements occurring during a sojournin state j . The rewards collected in each step can then be written as τ k = N (cid:88) j =0 τ k ( j )1 { J k = j } and ξ k = N (cid:88) j =0 ξ k ( j )1 { J k = j } . In the statements of our main theorems it will be useful to have a notationfor a vector of random variables with distributions for the time durations andspatial displacements that are associated with the states S = { , , . . . , N } : τ D = ( τ (0) , τ (1) , . . . , τ ( N )) and ξ D = ( ξ (0) , ξ (1) , . . . , ξ ( N )) . (14)So, for any step number k ∈ N , we have that the vector ( τ k (0) , τ k (1) , . . . , τ k ( N ))is equal in distribution to the vector τ and likewise for the spatial displace-ments. enewal reward perspective on linear switching diffusion systems 11 Example 1 (2-state advection-diffusion model of particle transport). Con-sider a 2-state advection-diffusion model for the dynamics of protein particles(such as mRNA) as illustrated in Ciocanel et al (2017, Figure 3A), with a freelydiffusing state and an active transport state. Assume that the times spent bythe particles in each state are drawn from an exponential distribution τ (0) ∼ Exp( β ) ,τ (1) ∼ Exp( β ) . Here β and β are the transition rates between states and the notation Exp( r )denotes an exponential distribution with parameter r (equal to the inverse ofthe mean). The spatial displacement in each state is given by: ξ (0) = (cid:112) Dτ (0) Z ,ξ (1) = vτ (1) , where D is the diffusion coefficient in the freely diffusing state, v is the speedin the active transport state, and Z are independent standard normal randomvariables. Example 2 (4-state reaction-diffusion-advection model of particle trans-port). More realistic representations of the dynamics of cellular protein con-centrations lead to considering the more complex 4-state model illustrated inCiocanel et al (2017, Figure 3B), where particles may diffuse, move in oppo-site directions, or be paused. The state durations are exponentially distributed,with the switching rates between dynamical states provided in Ciocanel et al(2017, Figure 3B), and the spatial displacements in each state are given by: ξ (0) = (cid:112) Dτ (0) Z ,ξ (1) = v + τ (1) ,ξ (2) = v − τ (2) ,ξ (3) = 0 , with v + the particle speed in the forward active transport state and v − theparticle speed in the backward active transport state. Example 3 (Cooperative models of cargo transport). Consider the coop-erative transport models proposed in Klumpp and Lipowsky (2005); Kunwarand Mogilner (2010), where processive motors move cargo in one directionalong a one-dimensional filament. These models assume a maximum number N of motor proteins, firmly bound to the cargo, that may act simultaneouslyin pulling the cargo in a specified direction (see Klumpp and Lipowsky (2005,Figure 1) for model visualization). The biophysical state (dynamic behavior)is defined by the number 0 ≤ n ≤ N of these motors that are bound to afilament and therefore actively contributing to transport. In a state with n motors attached to a filament, the cargo moves at a velocity v n , motors canunbind from the filaments with rate (cid:15) n or additional motors can bind to thefilaments with rate π n . The expressions for these transport model parametersare reproduced from Kunwar and Mogilner (2010), together with a nonlinearforce-velocity relation: v n ( F ) = v (cid:18) − (cid:18) FnF s (cid:19) w (cid:19) , (15) (cid:15) n ( F ) = n(cid:15)e F/ ( nF d ) , (16) π n = ( N − n ) π . (17)Here v is the load-free velocity of the motor, (cid:15) is the load-free unbinding rate,and π is the motor binding rate. F is the externally applied load force, F s is thestall force and F d is the force scale of detachment. The exponent w determinesthe nature of the force-velocity relation considered, with w = 1 correspondingto a linear relation, w < w > n (with 0 ≤ n ≤ N motors bound to the filaments) are therefore given by: τ ( n ) ∼ Exp (cid:0) r out ( n ) (cid:1) ,ξ ( n ) = v n ( F ) τ ( n ) , (18)where r out ( n ) = (cid:15) n ( F )+ π n is the transition rate out of the state with n motors(see Klumpp and Lipowsky (2005, Figure 1)). Example 4 (Tug-of-war models of cargo transport). Cargoes often movebidirectionally along filaments, driven by both plus and minus-directed motors.For example, kinesin moves cargo towards the plus end of microtubules whiledynein moves it towards the minus end. In M¨uller et al (2008, 2010), theauthors propose a model where a tug-of-war between motors drives cargo inopposite directions, with transport by several motors leading to an increasein the time the cargo remains bound to a microtubule and is pulled along aparticular direction. In these models, teams of maximum N + plus- and N − minus-end motors are bound to the cargo, and the biophysical state is givenby the pair of indices ( n + , n − ) with 0 ≤ n + ≤ N , 0 ≤ n − ≤ N indicatingthe number of plus and minus motors bound to the filament and therebycontributing actively to the transport (see M¨uller et al (2008, Figure 1) for enewal reward perspective on linear switching diffusion systems 13 model visualization). A key assumption for this model is that motors interactwhen bound to the filament since opposing motors generate load forces, andmotors moving in the same direction share the load. In addition, they assumethat motors move with the same velocity as the cargo in any state (M¨ulleret al, 2008, 2010). This model uses the following expressions for the transportparameters: v c ( n + , n − ) = n + F s + − n − F s − n + F s + /v f + + n − F s − /v b − , (19) (cid:15) + ( n + ) = n + (cid:15) e F c / ( n + F d + ) , (20) π + ( n + ) = ( N + − n + ) π . (21)Here indices + and − refer to the plus- and minus-end directed motors un-der consideration. The model parameters are as follows: F s is the stall force, F d is the force scale for detachment, (cid:15) is the load-free unbinding rate, π isthe motor binding rate, v f is the forward velocity of the motor (in its pre-ferred direction of motion), and v b is the slow backward velocity of the motorconsidered. Eq. 19 applies for the case when n + F s + > n − F s − (stronger plusmotors, M¨uller et al (2008)), and an equivalent expression with v f + replacedby v b + and v b − replaced by v f − holds for n + F s + ≤ n − F s − (stronger minusmotors). Equivalent expressions for the binding and unbinding rates hold forthe minus-end directed motors. In the case of stronger plus motors, the cargoforce F c when pulled by n + plus and n − minus motors is given by (M¨ulleret al, 2008): F c ( n + , n − ) = λn + F s + + (1 − λ ) n − F s − ,λ = 11 + n + F s + v b − n − F s − v f + , (22)with equivalent expressions for stronger minus motors as described above andin M¨uller et al (2008). The times and displacements accumulated at each timestep and in each state are defined as in Eq. 18 in Example 3. From standard renewal-reward and functional central limit theorem results,which we detailed in Section 2, we have related the computation of effectivevelocity and diffusivity via Eqs. 10 and 11 to analyzing the first and secondmoments and correlation of the spatial displacement and time spent in eachregeneration cycle. In this section, the main result is Proposition 1, whichprovides these statistics. We begin with Lemma 1, by recalling a standardrecursion formula for the moments of the reward accumulated until hitting adesignated absorbing state. We include the proof of this lemma for complete-ness and as an example of the moment generating function approach we use in Lemma 2. In Proposition 1 we address the calculation of total displacementand time duration during the regeneration cycles described in Section 2.Let 0 be the base state that marks the beginning of a new renewal cycle.We denote the set of remaining states as S \{ } , and define ˜ P as the N × N substochastic matrix containing the probabilities of transition among thesenon-base states only. Generally, we use the symbol ˜ to refer to a vector or amatrix whose components corresponding to the base state have been removed.Let R denote the total reward accumulated until the state process hits thebase state. Note that the value of R will depend on what the initial state of theprocess is. In our motor transport examples, R corresponds to the time ∆ ˜ T orthe displacement ∆ ˜ X accumulated after stepping away from the base state andbefore returning to the base state. Let ρ k denote the reward accumulated ateach time step, recalling that time increments are denoted τ k and displacementincrements ξ k in Section 2.2.By introducing random variables ρ k ( j ) for j ∈ S and k ∈ N that indicatethe reward received at step k if the particle is also in state j at that step, wecan use indicator variables for the state to express: ρ k = N (cid:88) j =1 ρ k ( j )1 { J k = j } and R = η (cid:88) k =1 ρ k = η (cid:88) k =1 N (cid:88) j =1 ρ k ( j )1 { J k = j } . (23)In the same way that we defined the distribution for the time durationsand spatial displacements through the random vectors τ and ξ in Eq. 14,we define the distribution of generic rewards through the vector of randomrewards associated to each state: ˜ ρ = ( ρ (1) , ρ (2) , . . . , ρ ( N )) . (24)The tilde notation is used here to be consistent with the connotation that tildeimplies the zero state is excluded. When we need component-wise multiplica-tion, we use the Hadamard power notation (see Eq. 5): ˜ ρ ◦ n = ( ρ n (1) , ρ n (2) , . . . , ρ n ( N )) . (25)We define the moment-generating functions of the reward collected untilthe state process hits the base state, and of the reward in state i , respectively,by the following vectors: φ ( s ) : φ i ( s ) := E ( e sR | J = i ) , and ψ ( s ) : ψ i ( s ) := E ( e sρ ( i ) ) . (26)Characteristic functions could alternatively be used to handle rewards whosehigher moments are not all finite; the results for the low order moments wecalculate would be the same. Note that here and in the following, we willtypically use index i to refer to states i ∈ S \{ } . In Lemma 1, J the statein the initial step of the process. We seek a general recursion relation for enewal reward perspective on linear switching diffusion systems 15 E ( R n | J = i ) and denote the corresponding vector of moments for all i ∈ S \{ } by E S \{ } ( R n ). The following result is a variation on similar recursionformulas for rewards accumulated in Markov chains (Hunter, 2008; Palacios,2009). Lemma 1
Let { J k } k ≥ be a time-homogeneous, positive recurrent Markovchain with a transition probability matrix P (over a finite state space S ) thathas zeroes for each of its diagonal entries. Let the reward variables R and ˜ ρ bedefined as in Eqs. 23 and 24, respectively. For n ∈ N , define the column vector E S \{ } ( R n ) := (cid:0) E ( R n | J = 1) , . . . , E ( R n | J = N ) (cid:1) . (27) Then this vector – the expected reward accumulated up to the first time thatthe state process { J k } hits the base state 0 – satisfies the recursion relation E S \{ } ( R ) = ( I − ˜ P ) − E ( ˜ ρ ); E S \{ } ( R n ) = ( I − ˜ P ) − (cid:32) E ( ˜ ρ ◦ n ) + n − (cid:88) m =1 (cid:18) nm (cid:19) diag( E ( ˜ ρ ◦ ( n − m ) )) ˜ P E S \{ } ( R m ) (cid:33) . (28) Here ˜ P is the substochastic matrix component of P excluding the base state ,and ˜ ρ ◦ n is the Hadamard n -th power vector defined in Eq. 25.Proof Let R , the reward accumulated until hitting the base state 0, be de-composed into the reward from the first and from subsequent steps as follows: R = ρ + ˇ R . We calculate the moment-generating function of R conditionedon the initial state J = i as follows: φ i ( s ) := E ( e sR | J = i )= (cid:88) j ∈ S E ( e sR | J = i, J = j ) P ij = (cid:88) j ∈ S E ( e sρ e s ˇ R | J = i, J = j ) P ij = E ( e sρ | J = i ) E ( e s ˇ R | J = 0) P i + (cid:88) j ∈ S \{ } E ( e s ˇ R | J = j ) P ij = E ( e sρ ( i ) ) P i + (cid:88) j ∈ S \{ } E ( e sR | J = j ) P ij = ψ i ( s ) P i + (cid:88) j ∈ S \{ } φ j ( s ) P ij , (29)where ψ i ( s ) is defined in Eq. 26. In the fourth line we used the Markov property,and in the fifth line we used the fact that( ˇ R | J = j ) ∼ ( R | J = j )(1 − δ j ) where δ ij is the Kronecker delta function. Defining f i ( s ) = ψ i ( s ) P i , i ∈ S \{ } ,G ( s ) = { G ( s, i, j ); i, j ∈ S \{ } : G ( s, i, j ) = ψ i ( s ) P ij } , then we can write Eq. 29 in matrix-vector form: φ ( s ) = f ( s ) + G ( s ) φ ( s ) . (30)Since the moments of the reward before hitting the base state can be calculatedusing E ( R n | J = i ) = ∂ n ∂s n φ i ( s ) | s =0 , we calculate the derivatives: ∂ n φ ( s ) ∂s n = ∂ n f ( s ) ∂s n + n (cid:88) m =0 (cid:18) nm (cid:19) ∂ n − m G ( s ) ∂s n − m ∂ m φ ( s ) ∂s m . For the first moment ( n = 1), each component yields: ∂φ i ( s ) ∂s = P i E ( ρ | J = i ) + (cid:88) j ∈ S \{ } P ij E ( ρ | J = i ) + (cid:88) j ∈ S \{ } P ij ∂φ j ( s ) ∂s = E ( ρ ( i )) + (cid:88) j ∈ S \{ } P ij ∂φ j ( s ) ∂s n . Evaluating at s = 0 for n = 1, we have E ( R | J = i ) = E ( ρ ( i )) + (cid:88) j ∈ S \{ } P ij E ( R | J n = j ) . Writing in vector form and solving for E S \{ } ( R ) yields the first part of Eq. 28.For higher-order moments ( n > ∂ n φ i ( s ) ∂s n = P i E ( ρ n | J = i ) + (cid:88) j ∈ S \{ } P ij E ( ρ n | J = i )+ n − (cid:88) m =1 (cid:18) nm (cid:19) E ( ρ n − m | J = i ) (cid:88) j ∈ S \{ } P ij ∂ m φ j ( s ) ∂s m + (cid:88) j ∈ S \{ } P ij ∂ n φ j ( s ) ∂s n = E ( ρ ( i ) n ) + (cid:88) j ∈ S \{ } P ij ∂ n φ j ( s ) ∂s n + n − (cid:88) m =1 (cid:18) nm (cid:19) E ( ρ ( i ) n − m ) (cid:88) j ∈ S \{ } P ij ∂ m φ j ( s ) ∂s m . Evaluating at s = 0 gives the recursion relation expressed in the second partof Eq. 28. enewal reward perspective on linear switching diffusion systems 17 Corollary 1
Let τ and ξ denote the vectors of state-dependent time durationand spatial displacements as defined in Eq. 14. Let ∆T and ∆X denote thetotal time elapsed and displacement accumulated by a state-switching particleup until its state process { J k } k ≥ returns to the base state 0 (see Eqs. 12).Moreover, recall the first-step decomposition ∆T = τ + ∆ ˜ T and ∆X = ξ + ∆ ˜ X (see Eqs. 13). Suppose that the state process { J k } k ≥ and its associatedtransition probability matrix P satisfy the assumptions of Lemma 1. Then E S \{ } ( ∆ ˜ T ) = ( I − ˜ P ) − E ( ˜ τ ) , (31) E S \{ } ( ∆ ˜ X ) = ( I − ˜ P ) − E ( ˜ ξ ) , (32) E S \{ } ( ∆ ˜ T ) = ( I − ˜ P ) − (cid:16) E ( ˜ τ ◦ ) + 2diag( E ( ˜ τ )) ˜ P E S \{ } ( ∆ ˜ T ) (cid:17) (33) E S \{ } ( ∆ ˜ X ) = ( I − ˜ P ) − (cid:16) E ( ˜ ξ ◦ ) + 2diag( E ( ˜ ξ )) ˜ P E S \{ } ( ∆ ˜ X ) (cid:17) , (34)where ˜ τ and ˜ ξ are the vectors of time durations and spatial displacementsexcluding the base state. Proof
These results follow directly from Lemma 1, with ∆ ˜ T and ∆ ˜ X respec-tively playing the role of the reward R . Lemma 2
Let τ , ξ , ∆ ˜ T , ∆ ˜ X , P , and { J k } k ≥ be defined as in Corollary 1.Then E S \{ } ( ∆ ˜ T ∆ ˜ X ) = ( I − ˜ P ) − (cid:16) E ( ˜ τ ◦ ˜ ξ ) + diag( E ( ˜ ξ )) ˜ P E S \{ } ( ∆ ˜ T )+ diag( E ( ˜ τ )) ˜ P E S \{ } ( ∆ ˜ X ) (cid:17) . (35) Proof
We use an argument similar to the re-arrangement of the moment-generating function in Eq. 30 in the proof of Lemma 1. Here we decomposethe time and displacement into the first step after the base state and the sub-sequent steps: ∆ ˜ T = τ + ˇ T and ∆ ˜ X = ξ + ˇ X . Since we are interested inthe cross-moment of the duration and displacement, we consider the followingmoment-generating function: φ i ( s, r ) = E ( e s∆ ˜ X e r∆ ˜ T | J = 0 , J = i )= (cid:88) j ∈ S E ( e s∆ ˜ X e r∆ ˜ T | J = 0 , J = i, J = j ) P ij = (cid:88) j ∈ S E ( e sξ e rτ e s ˇ X e r ˇ T | J = i, J = j ) P ij = E ( e sξ e rτ | J = i ) (cid:88) j ∈ S E ( e s ˇ X e r ˇ T | J = i, J = j ) P ij = E ( e sξ ( i ) e rτ ( i ) ) (cid:88) j ∈ S E ( e s ˇ X e r ˇ T | J = j ) P ij = ψ i ( s, r ) P i + ψ i ( s, r ) (cid:88) j ∈ S \{ } φ j ( s, r ) P ij , (36) where ψ i ( r, s ) = E ( e sξ ( i ) e rτ ( i ) ).For the calculation of the cross-term E S \{ } ( ∆ ˜ T ∆ ˜ X ), we note that ∂ φ i ∂r∂s | s = r =0 = E ( ∆ ˜ T ∆ ˜ X | J = i ) and calculate: ∂ φ i ∂r∂s = ∂ ∂r∂s ψ i ( r, s ) (cid:88) j ∈ S \{ } φ j ( r, s ) P ij + ψ i ( r, s ) P i = ∂∂r ∂ψ i ( r, s ) ∂s (cid:88) j ∈ S \{ } φ j ( r, s ) P ij + ψ i ( r, s ) (cid:88) j ∈ S \{ } ∂φ j ( r, s ) ∂s P ij + ∂ψ i ( r, s ) ∂s P i = ∂ ψ i ( r, s ) ∂s∂r (cid:88) j ∈ S \{ } φ j ( r, s ) P ij + ∂ψ i ( r, s ) ∂s (cid:88) j ∈ S \{ } ∂φ j ( r, s ) ∂r P ij + ∂ψ i ( r, s ) ∂r (cid:88) j ∈ S \{ } ∂φ j ( r, s ) ∂s P ij + ψ i ( r, s ) (cid:88) j ∈ S \{ } ∂ φ j ( r, s ) ∂s∂r P ij + ∂ ψ i ( r, s ) ∂s∂r P i . Evaluating the above at s = r = 0 yields: E ( ∆ ˜ T ∆ ˜ X | J = i ) = E ( τ ξ | J = i ) + E ( ξ | J = i ) (cid:88) j ∈ S \{ } P ij E ( ∆ ˜ T | J = j )+ E ( τ | J = i ) (cid:88) j ∈ S \{ } P ij E ( ∆ ˜ X | J = j )+ (cid:88) j ∈ S \{ } P ij E ( ∆ ˜ T ∆ ˜ X | J = j ) . Therefore, E S \{ } ( ∆ ˜ T ∆ ˜ X ) = E ( ˜ τ ◦ ˜ ξ ) + diag( E ( ˜ ξ )) ˜ P E S \{ } ( ∆ ˜ T )+ diag( E ( ˜ τ )) ˜ P E S \{ } ( ∆X ) + ˜ P E S \{ } ( ∆T ∆X ) , which yields equation 35. Remark 2
An alternative derivation of equation 35 would be to use a polar-ization argument for the expectation of the product: E S \{ } ( ∆ ˜ T ∆ ˜ X ) = 14 (cid:16) E S \{ } (( ∆ ˜ X + ∆ ˜ T ) ) − E S \{ } (( ∆ ˜ X − ∆ ˜ T ) ) (cid:17) . (37)In this approach, the moment-generating function depending on both cycletime ∆ ˜ T and cycle displacement ∆ ˜ X introduced in Eq. 36 is not required,since Lemma 1 can be directly applied to give explicit formulas for the secondmoments of the reward R = ∆ ˜ X + ∆ ˜ T and R = ∆ ˜ X − ∆ ˜ T .We proceed to Proposition 1, which provides the quantities necessary tocompute the effective velocity and diffusivity of the cargo dynamics usingclassical theory (see Eqs. 10 and 11 and the procedure in § enewal reward perspective on linear switching diffusion systems 19 Proposition 1 (First and second order statistics of rewards in a re-newal cycle)
Consider a regenerative cycle of a discrete-time time-homogeneous recur-rent Markov chain which takes its values in the discrete state space S = { , , , . . . , N } with probability transition matrix P with zero diagonal entries,starting at base state until its first return to base state . The associated time ∆T and spatial displacement ∆X are defined as in Eq. 12. The random vari-ables τ (0) and ξ (0) have the distributions of the time duration and spatialdisplacement that are accumulated in the base state, and p (1) is the vector oftransition probabilities from the base state in the first step of a cycle, i.e. thefirst row of P . The moments of the cycle time and displacement rewards arethen given by: E ( ∆T ) = E ( τ (0)) + p (1) · E S \{ } ( ∆ ˜ T ) , E ( ∆X ) = E ( ξ (0)) + p (1) · E S \{ } ( ∆ ˜ X ) , Var( ∆T ) = Var( τ (0)) + p (1) · E S \{ } ( ∆ ˜ T ) − ( p (1) · E S \{ } ( ∆ ˜ T )) , Var( ∆X ) = Var( ξ (0)) + p (1) · E S \{ } ( ∆ ˜ X ) − ( p (1) · E S \{ } ( ∆ ˜ X )) , Cov( ∆X, ∆T ) = Cov( τ (0) , ξ (0)) + p (1) · E S \{ } ( ∆ ˜ T ∆ ˜ X ) − ( p (1) · E S \{ } ( ∆ ˜ T ))( p (1) · E S \{ } ( ∆ ˜ X )) , (38) where the first, second, and cross-moments of the time ∆ ˜ T and the displace-ment ∆ ˜ X are given by Eqs. 31, 32, 33, 34, 35 in Corollary 1 and Lemma 2.Proof With state 0 as base state, we decompose the cycle time into the timespent in the base state τ = τ (0) and the time ∆ ˜ T spent from leaving the basestate until returning to the base state. Therefore, the total time in a cycle isgiven by ∆T = τ + ∆ ˜ T , and similarly the total spatial displacement in a cycleis ∆X = ξ + ∆ ˜ X . We apply the law of total expectation by conditioning onthe state J that the process visits after the base state: E ( ∆T ) = E ( E ( ∆T | J ))= (cid:88) i ∈ S \{ } E ( ∆T | J = i ) P i = (cid:88) i ∈ S \{ } E ( τ + ∆ ˜ T | J = i ) P i = E ( τ (0)) + (cid:88) i ∈ S \{ } E ( ∆ ˜ T | J = i ) P i = E ( τ (0)) + p (1) · E S \{ } ( ∆ ˜ T ) , (39)where as before S \{ } is the set of transient states and P i is the probabilityof switching from base state 0 to state i . A similar calculation applies to the first moment of the cycle reward E ( ∆X ). For the second moments, we use thelaw of total variance as follows:Var( ∆T ) = E (Var( ∆T | J )) + Var( E ( ∆T | J ))= E (Var( τ + ∆ ˜ T | J )) + Var( E ( τ + ∆ ˜ T | J ))= E (Var( τ (0)) + Var( ∆ ˜ T | J )) + Var( E ( τ (0)) + E ( ∆ ˜ T | J ))= Var( τ (0)) + E (Var( ∆ ˜ T | J )) + Var( E ( ∆ ˜ T | J ))= Var( τ (0)) + (cid:88) i ∈ S \{ } Var( ∆ ˜ T | J = i ) P i + (cid:88) i ∈ S \{ } ( E ( ∆ ˜ T | J = i )) P i − (cid:88) i ∈ S \{ } E ( ∆ ˜ T | J = i ) P i = Var( τ (0)) + (cid:88) i ∈ S \{ } E ( ∆ ˜ T | J = i ) P i − (cid:88) i ∈ S \{ } E ( ∆ ˜ T | J = i ) P i = Var( τ (0)) + p (1) · E S \{ } ( ∆ ˜ T ) − ( p (1) · E S \{ } ( ∆ ˜ T )) , (40)and similarly for Var( ∆X ). The covariance term can then be obtained via thepolarization formula from the formulas for the variances. Proposition 1 and the calculation procedure in § § P = (cid:18) (cid:19) . We take the diffusing state as the base state. The substochastic matrixof the probabilities of transition between the other states (Bhat and Miller,2002) is then simply the scalar ˜ P = 0 in this case, while the vector of transitionout of the base state is simply p (1) = [1] . enewal reward perspective on linear switching diffusion systems 21 The first and second moments of the cycle duration are given by Eqs 31and 33 with ( I − ˜ P ) − = 1. Similarly, the moments of the spatial displacementare given by Eqs. 32 and 34. In this model, we have that S \{ } = { } and˜ τ ◦ nk = τ nk (1) for the time reward and ξ ◦ nk = ξ nk (1) for the spatial displacementreward in the active transport state. In the 2-state system, these values aresimply scalars: E ( ∆ ˜ T ) = E ( τ (1)) = 1 /β , E ( ∆ ˜ X ) = E ( ξ (1)) = v/β , E ( ∆ ˜ T ∆ ˜ X ) = E ( τ (1) ξ (1)) = 2 v/β , E ( ∆ ˜ T ) = E ( τ (1) ) = 2 /β , E ( ∆ ˜ X ) = E ( ξ (1) ) = 2 v /β . The statistics of the cycle are therefore given by: E ( ∆T ) = E ( τ (0)) + E ( ∆ ˜ T ) p (1) (1) = 1 β + 1 β , E ( ∆X ) = E ( ξ (0)) + E ( ∆ ˜ X ) p (1) (1) = 0 + v/β = vβ , Var( ∆T ) = Var( τ (0)) + E ( ∆ ˜ T ) p (1) (1) − ( E ( ∆ ˜ T ) p (1) (1)) = 1 β + 1 β , Var( ∆X ) = Var( ξ (0)) + E ( ∆ ˜ X ) p (1) (1) − ( E ( ∆ ˜ X ) p (1) (1)) = 2 Dβ + v β , Cov( ∆T, ∆X ) = Cov( τ (0) , ξ (0)) + E ( ∆ ˜ T ∆ ˜ X ) p (1) (1) − (cid:16) E ( ∆ ˜ T ) p (1) (1) (cid:17) (cid:16) E ( ∆ ˜ X ) p (1) (1) (cid:17) = vβ . Eqs. 10 and 11 then provide expressions for the effective velocity and diffusivityof the particles as in Hughes et al (2011, 2012); Whitt (2002): v eff = E ( ∆X ) E ( ∆T ) = v β β + β ,D eff = 12 E ( ∆T ) ( v Var( ∆T ) + Var( ∆X ) − v eff cov( ∆T, ∆X ))= D β β + β + v β β ( β + β ) . Note that the effective velocity is given by the speed in the transportstate multiplied by the fraction of time the mRNA particles spend in themoving state. The effective diffusivity has a more complicated expression, butclearly shows the dependence of this quantity on each model parameter. Theseexpressions agree with the results of Eqs. 3, 4 as outlined in Ciocanel (2017). I − ˜ P (see Bhat and Miller (2002); Dobrow (2016)) to determine thefundamental matrix, the approach presented here is easily implemented in asoftware package such as Mathematica or Matlab for symbolic derivation ofthe effective transport properties for models with multiple states (see samplecode in the repository on GitHub (2019)).In Figure 2, we illustrate the good agreement of the results in Ciocanelet al (2017) with our calculation procedure in § N R = 500 stochastic realizations of the dynamics and foreach iteration, we run the process until a fixed large time T f = 5 × , whichalso keeps the computation feasible. We then estimate the effective velocityand diffusivity as follows: v eff ≈ ( (cid:80) N R i =1 X i ( T f )) /N R T f , (41) D eff ≈ (cid:16)(cid:80) N R i =1 ( X i ( T f ) − ( (cid:80) N R i =1 X i ( T f )) /N R ) (cid:17) / ( N R − T f , (42)where X i ( T f ) are the simulated final positions of the particle at time T f initeration i .The different parameter sets (labeled by index) in Figure 2 correspond tosimulations using parameter estimates based on FRAP mRNA data from dif-ferent frog oocytes in Gagnon et al (2013); Ciocanel et al (2017). The goodagreement of the theoretical and simulated effective velocity and diffusivityshows that the analytical approach proposed is a good alternative to poten-tially costly simulations of the stochastic process up to a large time. The framework presented here also extends to models of cargo particles drivenby changing numbers of motor proteins. The analytical calculation of trans- enewal reward perspective on linear switching diffusion systems 23
Parameter set index E ff e c t i v e v e l o c i t y ( m / s ) PDE analysisEstimate from simulationsRenewal rewards analysis -0.1-0.0500.050.10.15 E ff e c t i v e d i ff u s i v i t y ( / s ) m Parameter set index
PDE analysisEstimate from simulationsRenewal rewards analysis BA Fig. 2
Effective velocity (A) and effective diffusivity (B) of particles switching between dif-fusion, bidirectional transport, and stationary states as in Ciocanel et al (2017, Figure 3B)for different parameter sets. Blue triangles correspond to predictions based on the homog-enization or equivalent analysis (Ciocanel et al, 2017) of the corresponding PDEs (Eq. 1),filled red dots correspond to estimates from multiple simulated realizations of the Markovchain, and yellow circles correspond to predictions based on analysis of the correspondingrenewal process model combined with Proposition 1. port properties of cargo pulled by motors in the same or opposite directionscould replace or complement costly numerical simulations of individual cargotrajectories. In the following, we consider both models of cooperative cargotransport with identical motors (Klumpp and Lipowsky, 2005; Kunwar andMogilner, 2010) and tug-of-war models of bidirectional transport driven byidentical or different motors moving in opposite directions (M¨uller et al, 2008,2010).5.1 Cooperative models of cargo transportWe start by considering the cooperative transport models described in §
2, Ex-ample 3, and studied in Klumpp and Lipowsky (2005); Kunwar and Mogilner(2010), with processive motors that move along a one-dimensional microtubuleand transport cargo in only one direction. The dynamics is described by theforce-driven velocities v n , unbinding rates (cid:15) n and binding rates π n in each statewith n motors bound to the microtubule and moving the cargo (see Eqs. 15,16, and 17). In this section, we use the kinetic parameters for conventionalkinesin-1 provided in Klumpp and Lipowsky (2005).Our calculation of the effective velocity of cargo agrees with the derivationin Klumpp and Lipowsky (2005), which uses the stationary solution of themaster equation for probabilities of the cargo being in each state (i.e. carriedby n motors). We note that there are two notions of effective velocity (and dif-fusivity) that can be used in studying this model: one is to calculate the effec-tive velocity of the cargo averaged over the bound states only (the asymptoticvelocity without detachment along a theoretical infinite length microtubule)(Klumpp and Lipowsky, 2005; Kunwar and Mogilner, 2010), and the secondis to calculate the overall effective velocity that also accounts for periods of detachment from microtubules. For the N = 2 motors model, Klumpp andLipowsky (2005) and Kunwar and Mogilner (2010) report the average velocityfor bound cargo (first notion): v eff = v π (cid:15) π (cid:15) + π π + v π π π (cid:15) + π π . (43)Since we are interested in the overall effective velocity of the particles in thecontext of their full dynamics, we include the state where no motors are boundto the filament in our calculation, so that the effective velocity with respectto the overall dynamics is given by: v eff = v π (cid:15) (cid:15) (cid:15) + π (cid:15) + π π + v π π (cid:15) (cid:15) + π (cid:15) + π π . (44)Using the calculation of the overall effective velocity in (44), we predicta similar dependence of the effective velocity under a range of force loads asin Klumpp and Lipowsky (2005) using the formula (43). The dashed curvesin Figure 3A,C agree with the behavior of sub- and super-linear motors underdifferent load forces as reported in Kunwar and Mogilner (2010, Figure 2C-D),including the fact that sub-linear motors have lower effective velocities for anychoice of the load force and for all maximum motor numbers N considered(A, w = 0 . w = 2).The insight from our method lies in the prediction of the effective diffu-sivity as a function of load for each type of motor. Figure 3B,D show thatthe N = 1 motor transport case has a large effective diffusivity under no loadbecause of the switching between the paused and moving states. As the forceload increases to stall F s , the velocity of the single motor state decreases to0: v ( F ) = v (1 − F/F s ). Therefore, the active transport state switches to astationary state at F = F s = 6 pN, leading to decreased effective diffusivityas the cargo switches between dynamic states with similar behaviors.For N = 2 and N = 3, the calculation of the effective diffusivity allows usto re-visit the cooperative transport models for a large range of load forces andobserve a new phenomenon in the classical models of Klumpp and Lipowsky(2005); Kunwar and Mogilner (2010). The broader sweep of the load forceparameter in Figure 3B,D shows a non-monotonic dependence of the effectivediffusivity on load force for all types of motors considered (linear, sublinear,and superlinear), with an increase in effective diffusivity of cargo at low loadforces and a decrease at large load forces. While it is not immediately clearwhat leads to this phenomenon, we conjecture that this observation may bea result of the balance between two competing effects: on the one hand, asthe load increases, there is more detachment of motors (see (16)) and thusmore frequent switches between transport and stationary states, leading to anincrease in effective diffusivity; on the other hand, the increase in load forceleads to a decrease in the speeds of the motor-driven cargo states (see (15)),and thus a decrease in effective diffusivity. enewal reward perspective on linear switching diffusion systems 25 N=1, w=1N=2, w=1N=3, w=1N=1, w=2N=2, w=2N=3, w=2 E ff e c t i v e v e l o c i t y ( ) µ m / s Load force (pN)
N=1, w=1N=2, w=1N=3, w=1N=1, w=0.5N=2, w=0.5N=3, w=0.5 E ff e c t i v e v e l o c i t y ( ) µ m / s Load force (pN)
N=1, w=1N=2, w=1N=3, w=1N=1, w=2N=2, w=2N=3, w=2
Load force (pN) E ff e c t i v e d i ff u s i v i t y ( ) µ m / s N=1, w=1N=2, w=1N=3, w=1N=1, w=0.5N=2, w=0.5N=3, w=0.5
Load force (pN) E ff e c t i v e d i ff u s i v i t y ( ) µ m / s BAC D
Fig. 3
Effective velocity (A,C) and effective diffusivity (B,D) of cargo driven by a maximalnumber N of forward motor proteins as a function of the load force under various velocity-force exponents w (Eq. 15). The stall force used for kinesin is F s = 6 pN. Solid linescorrespond to motors with a linear force-velocity relation and dashed lines correspond tosub-linear motors with a convex-up force-velocity relation (top row), and respectively tosuper-linear motors with a concave force-velocity relation (bottom row). §
2, we consider the case where plus- and minus-directedmotors can drive cargo bidirectionally along filaments. The cargo velocities v c ( n + , n − ), unbinding rates (cid:15) + / − ( n + , n − ) and binding rates π + / − ( n + / − ) de-pend on the number of plus motors n + and minus motors n − at each state. Identical plus and minus motors.
With kinesin parameters drawn fromM¨uller et al (2008, Table 1), we first calculate the transport properties ofcargo in these models for identical plus and minus motors in equal numbers( N + = N − ). We vary the stall force of the kinesin motor to determine if thetheoretical effective velocity and diffusivity capture the differences obtained inthe numerical simulation studies in M¨uller et al (2008, 2010) for weak motors(small stall to detachment force ratio f = F s /F d ) and strong motors (large f ).As expected, the effective velocity in this symmetric case of identical motorsis zero for all stall forces (see Figure 4A). The predicted effective diffusivity inFigure 4B shows that for weak motors, the effective diffusivity is small, and dif- E ff e c t i v e d i ff u s i v i t y ( ) µ m / s N=1N=2N=3
Stall force (pN)
N=1N=2N=3 E ff e c t i v e d i ff u s i v i t y ( ) µ m / s Stall force (pN)
N=1N=2N=3 E ff e c t i v e v e l o c i t y ( ) µ m / s Stall force (pN)
N=1N=2N=3 E ff e c t i v e v e l o c i t y ( ) µ m / s Stall force (pN)
BAC D
Fig. 4
Effective velocity (A,C) and effective diffusivity (B,D) of cargo driven by maximum N forward and maximum N backward motor proteins as a function of kinesin-1 stall force F s ; the detachment force is F d = 3 pN. Panels (A,B) correspond to identical forward andbackward motors with kinetic parameters for kinesin-1, and panels (C,D) correspond tokinesin-1 forward motors and conventional dynein backward motors. ferent maximum numbers of motors do not lead to significant differences. Thisis similar to the results in M¨uller et al (2008), where the simulated cargo trajec-tories show small fluctuations and the probability distribution for the velocityhas a single maximum peak corresponding to approximately equal numbers ofplus and minus motors attached. However, for strong motors with a larger stallto detachment force ratio, the effective diffusivity increases considerably forall models. This is consistent with the observation in M¨uller et al (2008) thatstrong motors lead to cascades of unbinding of minus motors until only plusmotors stay bound (and vice versa), so that the spread of the cargo positionis predicted to be larger. The larger motor numbers lead to a more significantincrease in effective diffusivity as observed in M¨uller et al (2010), where thesimulated diffusion coefficient grows exponentially with motor numbers andtherefore leads to a more productive search of target destinations throughoutthe domain (M¨uller et al, 2010).It is worth noting that the method we develop in § § enewal reward perspective on linear switching diffusion systems 27 of the cargo does not change, however the effective diffusivity is consistentlylarger than in the case where the unbound cargo is fully stationary (resultsnot shown). Distinct plus and minus motors.
When considering dynein as the minus-end directed motor in the bidirectional transport model, we use the kineticparameters estimated to fit
Drosophila lipid droplet transport in M¨uller et al(2008, Table 1). Figure 4C shows that the cargo is predicted to move in theforward (kinesin-driven) direction with a positive effective velocity. We againobserve increased transport efficiency for larger numbers of motors. With in-creasing stall force, the velocity of individual runs in each state increases andtherefore the effective velocity increases and then plateaus. This asymmetricmotor case also results in effective diffusivity that decreases past a small stallforce and then stabilizes (see Figure 4D). Since the kinesin motor dominatesthe dynamics, there are fewer excursions backwards than in the case of iden-tical motors, so that the effective diffusivity is an order of magnitude smaller.Larger teams of motors regularize the dynamics and display decreased effectivediffusivity.5.3 Reattachment in models of cargo transportIn vitro experiments have suggested that binding rates of molecular motorsat specific locations may be regulated by the concentration of the same oropposite-directed motors (Hancock, 2014), as well as by the availability ofmicrotubule filaments. To test for the impact of reattachment kinetics in thestandard transport models of M¨uller et al (2008, 2010), we modify the bindingrate in (21) to account for a higher likelihood of reattachment when a motor(of either type) is already attached to the microtubule: π + ( n + , n − ) = (cid:40) N + π , if n + + n − = 0 , ( N + − n + ) ρπ , else . (45)Here ρ > π − ( n + , n − ). ρ = 1 corresponds tothe binding kinetics in the previous sections, and ρ > ∆X , panels C,D) for values of ρ ranging from1 to 50, in the context of models labeled ( N , N ) with transport driven bymaximum N forward (kinesin-1) motors and N backward (dynein) motors.Here we report the overall effective velocity of the cargo according to thesecond definition in § E ( ∆X ) in Eqs. 38. Notethe base state of complete detachment makes no contribution to the meandisplacement. (reattachment factor) E x pe c t ed r un l eng t h ( m ) (1,0)(1,1)(2,0)(2,1)(2,2) C (reattachment factor) E x pe c t ed r un l eng t h ( m ) (1,0)(1,1)(2,0)(2,1)(2,2) D (reattachment factor) E ff e c t i v e v e l o c i t y ( m / s ) (1,0)(1,1)(2,0)(2,1)(2,2) Kin1-DDBB (reattachment factor) E ff e c t i v e v e l o c i t y ( m / s ) (1,0)(1,1)(2,0)(2,1)(2,2) A Kin1-Dynein
Fig. 5
Effective velocity (A,B) and expected run length (C,D) of cargo driven by maximum N forward (kinesin-1) and maximum N backward (dynein) motor proteins ( N , N ) as afunction of the reattachment factor ρ . Panels (A,C) use conventional dynein motor kinetics asin M¨uller et al (2008, 2010) while panels (B,D) use dynein-dynactin-BicD2 (DDB) complexparameters as in Ohashi et al (2019). The run lengths in panels (B) and (D) are plotted ona log-log scale to allow for visualization of differences between the models considered. Thedefinitions of effective velocity and run length used are provided in the text. Classically in tug-of-war modeling, dynein has been viewed as a “weakerpartner” than kinesin-family motors. In this parameter regime (M¨uller et al,2008, 2010), dynein has both a smaller stall force and smaller critical detach-ment force than kinesin-1. As a result, when equal numbers of kinesin-1 anddynein are simultaneously attached, kinesin-1 dominates transport. However,it has recently been shown that it might not be realistic to consider dynein inthe absence of its helper proteins, particularly dynactin and BicD2. Together,these form a complex referred to as DDB, and the associated parameter val-ues (Ohashi et al, 2019) are much more “competitive” with kinesin-1 in atug-of-war scenario.In Figure 5, we display the effective velocity and expected run length ofkinesin-1 vs dynein (panels A,C), and kinesin-1 versus DDB (panels B,D) dy-namics. In Figure 5A, the effective velocity of cargo driven by teams of motorsapproaches the effective speed predicted for kinesin-only motor teams (mod-els (1 ,
0) and (2 , ρ , but then decreases as ρ becomeslarger for conventional dynein motility. As observed in recent studies, acti- enewal reward perspective on linear switching diffusion systems 29 Estimate from simulationsFirst passage approximationRenewal rewards asymptotics Load sensitivity
Log o f e ff e c t i v e d i ff u s i v i t y ( / s ) m Fig. 6
Comparison of effective diffusivity estimates for parallel microtubules driven bidirec-tionally by kinesin motors as a function of the scaled load sensitivity γ : stochastic simulationsfrom Allard et al (2019) marked with blue stars, first passage time approximation in Allardet al (2019) marked with red circles, and renewal rewards calculation marked with yellowtriangles. Following Allard et al (2019), we only allow states with i kinesin motors movingforward and K − i kinesin motors moving backward, with 0 ≤ i ≤ K and K = 35. Thevertical axis is plotted on a log scale to allow for visualization of differences between effectivediffusivity estimated from simulations and analytical approximations. vated dynein competes more efficiently with kinesin and therefore the teamsof opposite-directed motors are consistently slower than teams consisting ofonly the forward kinesin motor protein in Figure 5B (Ohashi et al, 2019). Theexpected run lengths in Figure 5C,D illustrate that teams of multiple mo-tors are characterized by significantly increased processivity on microtubulesas the reattachment factor becomes larger. When considering conventionaldynein, the difference in processive cargo motion between the cooperative andtug-of-war models is only observed at large values of the reattachment con-stant ρ ( >
10, see Figure 5C). This is due to the fact that the backward motor(conventional dynein) in the M¨uller et al (2008) model is weak with a smalldetachment force, so that overcoming the large dynein unbinding rate requireslarge values of the reattachment factor. On the other hand, activated dyneinin the DDB complex is a more equal competitor to kinesin, with predictions ofthe expected run length in Figure 5D confirming the experimental observationsof larger unloaded run lengths in Ohashi et al (2019).5.4 Microtubule sliding modelAs a final example of the applicability of our method, we consider a recentinvestigation into microtubule motility and sliding by Allard et al (2019). Theauthors consider a continuous-time Markov chain model of the interactionof two parallel microtubules, cross-linked and moved by multiple identicalkinesin motors. Depending on which microtubule the motor heads are attachedto, they push the microtubule pair apart in one of two directions (one ofwhich is arbitrarily assigned to be “positive”). The model assumes that motorattachment to microtubules occurs quickly relative to detachment, allowing a reduced number of dynamic states with microtubules driven by i motorspushing in the positive direction and K − i motors pushing in the negativedirection (where K is the maximal number of motors that fit the overlap regionbetween the two parallel microtubules). The detachment rates are thereforegiven by κ + i = ( K − i ) κ exp (cid:18) γ iK (cid:19) ,κ − i = iκ exp (cid:18) γ K − iK (cid:19) , (46)where κ + i is the rate at which a motor pulling in the negative direction isreplaced by one pulling in the positive direction, κ − i is the is the rate at whicha motor pulling in the positive direction is replaced by one moving in thenegative direction, κ is a force-free transition rate, and γ is a dimensionlessload sensitivity defined as twice the stall force divided by the detachment forcescale (Allard et al, 2019). The relative velocity of the parallel microtubules ineach state is given by ∆v i = V m i − KK , (47)where V m is the speed of a single motor (Allard et al, 2019).A main point of this study is that parallel microtubules may slide bidi-rectionally with respect to each other, with a zero mean velocity due to sym-metry, thus the long-term microtubule transport is characterized by diffusivebehavior. The effective diffusivity of the microtubule pair driven by a totalof 35 motors is measured in Allard et al (2019) through fitting the slope ofthe mean squared displacement in stochastic simulations at long time (starsin Figure 6) and compared to a theoretical approximation in terms of a firstpassage time problem (open circles in Figure 6).Our method based on renewal rewards theory yields predictions of the effec-tive diffusivity that are closer to the estimates derived from large-time stochas-tic simulations (marked with triangles in Figure 6). To make the relationshipbetween effective diffusivity and load sensitivity more clear, we illustrate re-sults for many intermediary values. Our proposed analytical framework alsofacilitates the possibility of subsequent systematic asymptotic approximationsto study dependence on underlying biophysical parameters. In this work, we consider examples from the intracellular transport literaturewhere particles undergo switching dynamics. In particular, we are interestedin determining the effective velocity and diffusivity as well as the expected runlength of these particles as they switch between biophysical behaviors such asdiffusion, active transport, and stationary states. We propose a method that enewal reward perspective on linear switching diffusion systems 31 is based on defining the underlying Markov chain of state switches and theindependent cycles of the dynamics marked by returns to a chosen base state.Emphasizing the cyclic structure of the behavior allows us to treat the timedurations and spatial displacements of particles in these regenerative cycles asthe cycle durations and rewards in a renewal-reward process. Through calcu-lation of the statistics of cycle time and displacement, this robust frameworkprovides a rigorous means to study how the dynamics of switching systemsdepends on model parameters.This approach applies for example to canonical tug-of-war models describ-ing the transport of cargo by teams of molecular motor proteins. Previousinvestigations of the effective transport of cargo in these multi-state modelshave considered individual trajectories of the dynamics, computed using MonteCarlo simulations with the Gillespie algorithm (M¨uller et al, 2008; Kunwarand Mogilner, 2010; M¨uller et al, 2010). These studies determine the effec-tive velocity of the particles analytically by calculating the distribution of thenumber of bound motors from the stationary solution of the master equation(Klumpp and Lipowsky, 2005). However, determining the effective diffusivityin these studies relied on numerical simulations. Our method proposes a fasterand explicit investigation of the impact of model parameters on the effectivediffusivity. For instance, Figure 4 (top right) captures the different behaviorof identical motor teams involved in tug-of-war dynamics when the ratio ofstall to detachment force is small (weak motors with small effective diffusivity)versus large (strong motors with increasing effective diffusivity). This obser-vation is consistent with simulations in M¨uller et al (2008), where the largeforce ratios correspond to a dynamic instability where only one motor type isprimarily bound at the end of an unbinding cascade (M¨uller et al, 2008, 2010).Multiple experiments summarized in Hancock (2014) have shown that in-hibition of one motor type reduces transport in both directions in severalsystems, suggesting a “paradox of co-dependence” in bidirectional cargo trans-port. Several mechanisms accounting for this paradox were proposed, includingthe microtubule tethering mechanism recently explored in Smith and McKin-ley (2018). The hypothesis for this mechanism is that motors switch betweendirected active transport and a weak binding or diffusive state. The recentexperimental study in Feng et al (2018) suggests that teams of kinesin-1 mo-tors coordinate transport using help from the dynamic tethering of kinesin-2 motors. This work shows that when kinesin-1 motors detach, tethering ofkinesin-2 to the microtubule ensures that cargo stays near the filament toallow for subsequent reattachment (Feng et al, 2018). Our approach allowsus to assess the dependence of the dynamics on a potentially increased reat-tachment rate for cargo that is already bound to the filament by at least onemotor (Figure 5). Implementing this change in the standard binding models inM¨uller et al (2008); Kunwar and Mogilner (2010); M¨uller et al (2010) for bothkinesin-1/dynein and kinesin-1/DDB dynamics shows a decrease in overall ef-fective velocity, but very large increases in potential run length. This couldbe consistent with the paradox in that experimentalists would observe morekinesin-directed activity when the reattachment rate is sufficiently high.
We have made Matlab and Mathematica sample code available for thecalculation of effective velocity, diffusivity, and run lengths in a cooperativemodel of one-directional transport (as discussed in 5.1) and a tug-of-war modelof bidirectional transport (as discussed in 5.2) (GitHub, 2019). The code forthese examples can be readily adapted to allow for a general probability transi-tion matrix for the state dynamics, together with the probability distributionsfor the times and displacement in each state, to extend to other models of theprocessive movement of molecular motors and cargo transport. As the theoryof how motors coordinate to transport cargo continues to develop at a rapidpace, the analysis developed here will provide a tool for new models account-ing for tethered and weakly binding states with stochastic transitions whoserates do not depend on spatial position.
References
Allard J, Doumic M, Mogilner A, Oelz D (2019) Bidirectional sliding of twoparallel microtubules generated by multiple identical motors. Journal ofMathematical Biology pp 1–24Bhat UN, Miller GK (2002) Elements of applied stochastic processes, vol 3.Wiley-Interscience Hoboken, NJBressloff PC, Newby JM (2011) Quasi-steady-state analysis of two-dimensionalrandom intermittent search processes. Physical Review E 83(6):061139Bressloff PC, Newby JM (2013) Stochastic models of intracellular transport.Reviews of Modern Physics 85(1):135Bressloff PC, Xu B (2015) Stochastic active-transport model of cell polariza-tion. SIAM Journal on Applied Mathematics 75(2):652–678Brooks EA (1999) Probabilistic methods for a linear reaction-hyperbolic sys-tem with constant coefficients. Annals of Applied Probability pp 719–731Ciocanel MV (2017) Modeling intracellular transport during messenger RNAlocalization in
Xenopus oocytes. PhD thesis, Brown UniversityCiocanel MV, Sandstede B, Jeschonek SP, Mowry KL (2018) Modelingmicrotubule-based transport and anchoring of mRNA. SIAM Journal onApplied Dynamical Systems 17(4):2855–2881Ciocanel V, Kreiling JA, Gagnon JA, Mowry KL, Sandstede B (2017) Analysisof active transport by fluorescence recovery after photobleaching. Biophys-ical Journal 112(8):1714–1725Cioranescu D, Donato P (1999) An Introduction to Homogenization. OxfordUniversity Press, New YorkCox DR (1962) Renewal theory. MethuenCox SM, Matthews PC (2002) Exponential time differencing for stiff systems.Journal of Computational Physics 176(2):430–455Dobrow RP (2016) Introduction to stochastic processes with R. John Wiley& Sons enewal reward perspective on linear switching diffusion systems 33
Feng Q, Mickolajczyk KJ, Chen GY, Hancock WO (2018) Motor reattach-ment kinetics play a dominant role in multimotor-driven cargo transport.Biophysical Journal 114(2):400–409Gagnon JA, Kreiling JA, Powrie EA, Wood TR, Mowry KL (2013) Directionaltransport is mediated by a dynein-dependent step in an RNA localizationpathway. PLOS Biol 11(4):e1001551GitHub (2019) Sample Matlab and Mathematica code for effective ve-locity and diffusivity calculation. https://github.com/scottmckinley/stochastics-lab/tree/master/effective-transport
Hancock WO (2014) Bidirectional cargo transport: moving beyond tug of war.Nature Reviews Molecular Cell Biology 15(9):615Hughes J, Hancock WO, Fricks J (2011) A matrix computational approach tokinesin neck linker extension. Journal of Theoretical Biology 269(1):181–194Hughes J, Hancock WO, Fricks J (2012) Kinesins with extended neck linkers:a chemomechanical model for variable-length stepping. Bulletin of Mathe-matical Biology 74(5):1066–1097Hunter JJ (2008) Variances of first passage times in a markov chain withapplications to mixing times. Linear Algebra and its Applications 429(5-6):1135–1162Jung P, Brown A (2009) Modeling the slowing of neurofilament transport alongthe mouse sciatic nerve. Physical Biology 6(4):046002Klumpp S, Lipowsky R (2005) Cooperative cargo transport by several molec-ular motors. Proceedings of the National Academy of Sciences of the UnitedStates of America 102(48):17284–17289Kramer PR, Latorre JC, Khan AA (2010) Two coarse-graining studies ofstochastic models in molecular biology. Communications in MathematicalSciences 8(2):481–517Krishnan A, Epureanu BI (2011) Renewal-reward process formulation of motorprotein dynamics. Bulletin of Mathematical Biology 73(10):2452–2482Kunwar A, Mogilner A (2010) Robust transport by multiple motors with non-linear force–velocity relations and stochastic load sharing. Physical Biology7(1):016012Kunwar A, Tripathy SK, Xu J, Mattson MK, Anand P, Sigua R, VershininM, McKenney RJ, Clare CY, Mogilner A, et al (2011) Mechanical stochas-tic tug-of-war models cannot explain bidirectional lipid-droplet transport.Proceedings of the National Academy of Sciences 108(47):18960–18965Lawler GF (1995) Introduction to stochastic processes. Chapman & Hall, NewYorkLi Y, Brown A, Jung P (2014) Deciphering the axonal transport kinetics ofneurofilaments using the fluorescence photo-activation pulse-escape method.BMC Neuroscience 15(Suppl 1):P132McKinley SA, Athreya A, Fricks J, Kramer PR (2012) Asymptotic analysis ofmicrotubule-based transport by multiple identical molecular motors. Journalof Theoretical Biology 305:54–69Miles CE, Lawley SD, Keener JP (2018) Analysis of nonprocessive molecularmotor transport using renewal reward theory. SIAM Journal on Applied
Mathematics 78(5):2511–2532M¨uller MJ, Klumpp S, Lipowsky R (2008) Tug-of-war as a cooperative mech-anism for bidirectional cargo transport by molecular motors. Proceedings ofthe National Academy of Sciences 105(12):4609–4614M¨uller MJ, Klumpp S, Lipowsky R (2010) Bidirectional transport by molecu-lar motors: enhanced processivity and response to external forces. Biophys-ical Journal 98(11):2610–2618Newby J, Bressloff PC (2010a) Random intermittent search and the tug-of-warmodel of motor-driven transport. Journal of Statistical Mechanics: Theoryand Experiment 2010(04):P04014Newby JM, Bressloff PC (2010b) Quasi-steady state reduction of molecularmotor-based models of directed intermittent search. Bulletin of Mathemat-ical Biology 72(7):1840–1866Ohashi KG, Han L, Mentley B, Wang J, Fricks J, Hancock WO (2019) Load-dependent detachment kinetics plays a key role in bidirectional cargo trans-port by kinesin and dynein. Traffic 20(4):284–294Palacios JL (2009) On the moments of hitting times for random walks on trees.Journal of Probability and Statistics 2009Pavliotis G, Stuart A (2008) Multiscale methods: averaging and homogeniza-tion. Springer Science & Business MediaPavliotis GA (2005) A multiscale approach to Brownian motors. Physics Let-ters A 344(5):331–345Popovic L, McKinley SA, Reed MC (2011) A stochastic compartmentalmodel for fast axonal transport. SIAM Journal on Applied Mathematics71(4):1531–1556Reed MC, Venakides S, Blum JJ (1990) Approximate traveling waves in lin-ear reaction-hyperbolic equations. SIAM Journal on Applied Mathematics50(1):167–180Reimann P, Van den Broeck C, Linke H, H¨anggi P, Rubi J, P´erez-Madrid A(2001) Giant acceleration of free diffusion by use of tilted periodic potentials.Physical review letters 87(1):010602Reimann P, Van den Broeck C, Linke H, H¨anggi P, Rubi J, P´erez-MadridA (2002) Diffusion in tilted periodic potentials: Enhancement, universality,and scaling. Physical Review E 65(3):031104Resnick S (1992) Adventures in stochastic processes. Birkh¨auser Boston Inc.,Boston, MASerfozo R (2009) Basics of applied stochastic processes. Springer Science &Business MediaShtylla B, Keener JP (2015) Mathematical modeling of bacterial track-alteringmotors: Track cleaving through burnt-bridge ratchets. Physical Review E91(4):042711Smith JD, McKinley SA (2018) Assessing the impact of electrostatic drag onprocessive molecular motor transport. Bulletin of Mathematical Biology pp1–36Trong PK, Doerflinger H, Dunkel J, St Johnston D, Goldstein RE (2015)Cortical microtubule nucleation can organise the cytoskeleton of
Drosophila enewal reward perspective on linear switching diffusion systems 35enewal reward perspective on linear switching diffusion systems 35