
==== Front
R Soc Open Sci
R Soc Open Sci
RSOS
royopensci
Royal Society Open Science
2054-5703
The Royal Society

rsos231588
10.1098/rsos.231588
1008659119Mathematics
Research Articles
A mathematical framework for analysing particle flow in a network with multiple pools
A mathematical framework for analysing particle flow in a network with multiple pools
Jain Aditi 1 Conceptualization Formal analysis Investigation Methodology Software Validation Writing – original draft 2018maz0007@iitrpr.ac.in

https://orcid.org/0000-0001-6671-6747
Gupta Arvind Kumar 1 Conceptualization Funding acquisition Methodology Project administration Supervision Writing – review and editing akgupta@iitrpr.ac.in

1 Department of Mathematics, Indian Institute of Technology Ropar , Rupnagar, Punjab 140001, India
5 2024
08 5 2024 May 8, 2024
08 5 2024 May 8, 2024
11 5 23158820 10 2023 October 20, 2023
30 1 2024 January 30, 2024
22 2 2024 February 22, 2024
© 2024 The Authors.
2024
https://creativecommons.org/licenses/by/4.0/ Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

In many real-world systems, the entry rate of particles into a lane is affected by the occupancy of nearby pools. For instance, in biological networks, the concentration of molecules on the side of a membrane affects the entry of particles through the membrane. To understand the behaviour of such networks, we develop a network model of ribosome flow models (RFMs) having multiple pools where each RFM captures the dynamics of particle flow in a lane and competes for the finite resources present at the nearby pool. We study a ribosome flow model network with two pools (RFMNTP) and show that the network always admits a steady state. We then analyse the behaviour of the RFMNTP with respect to modifying the transition rate through a theoretical framework. Simulations of the RFMNTP demonstrate a counterintuitive result. For example, increasing any of the transition rates in the presence of a slow site in an RFM can increase the output rate of some RFMs and decrease the output rate of the other RFMs simultaneously. This suggests that the role of local sharing of particles incorporated is non-trivial. Finally, we illustrate how these results can provide insights into studying a network with multiple pools.

network with multiple pools
; ribosome flow model
; cooperative dynamical system
; first integral
; entrainment
==== Body
pmc1. Introduction

Movement is a vital aspect of life. Various cellular and physical processes involve the movement of particles along tracks. These processes generally take place in parallel, and they compete for the available limited resources. For example, during gene expression, all DNA (mRNA) molecules simultaneously compete for the limited amount of RNAPs (ribosomes) [1–3]. This competition generates a network in which an indirect coupling induces interactions among the lanes even in the absence of explicit connections [4,5]. Hence, it is of considerable interest to analyse such networks in the presence of these interactions and also to design several resource-sharing synthetic gene networks [6,7]. In physical systems such as vehicular flow, the number of vehicles moving along the roads is finite. The entry rate of the vehicles along the road is affected owing to the queue of vehicles waiting to enter a road, where each vehicle competes with other vehicles for limited space on the road [8]. This requires modelling the complex road network to comprehend the traffic flow, thereby reducing travel time and preventing traffic deadlocks [9–11].

Several computational and mathematical models have been developed to study resource-sharing networks [12–14]. One such model includes a set of totally asymmetric simple exclusion processes (TASEPs) that are interconnected to each other via a pool of free particles [15–17]. In TASEP, a chain of sites models the one-dimensional track and particles hop unidirectionally, with some probability, from one site to the next site only if the target site is empty [18,19]. The model encapsulates the simple exclusion concept through the fact that particles have volume and are unable to pass a moving particle ahead. The TASEP and its networks have been used to model and analyse various natural and artificial systems, including mRNA translation, vehicular traffic flow and more [10,20,21]. Regardless of its simple description, rigorous analysis of networks of TASEPs is complex, exact solutions exist for simplified cases and most non-homogeneous cases are studied via numerical methods or extensive Monte Carlo simulations [22–24]. Hence, understanding the effect of parameters on the dynamics of a system through TASEP-based models has proved challenging. To allow analytical treatment, network models based on ribosome flow models (RFMs) were developed [25–27].

The RFM is a deterministic model that is derived via the mean-field approximation of TASEP [28,29]. It is a nonlinear model, and yet, it is very amenable to mathematical analysis using methods from systems and control theory, in particular, contraction theory [30–32], monotone dynamical systems [33,34], theory of continued fractions [35], etc. The TASEP-based models use mean-field approximations in the thermodynamic limit, and hence, analysis of the TASEP becomes accurate only when the number of sites n goes to infinity, while the analysis of the RFM provides results that are valid for every n . The theory of TASEP focuses on domain wall theory, phase transitions, etc. [36–38], which is different from the theoretical approach used to analyse RFM. The RFM and its generalizations have been developed to understand the complex dynamics of a variety of transport phenomena on a single lane, including translation [35,39–44], transcription [45], motor protein traffic [46,47], phosphorelay [48] and more. There are lattice hydrodynamic models that use ordinary differential equations to model the flow of vehicles along the lanes [49–51]. Therefore, the framework of RFM that also describes the flow of particles can serve as the basis for understanding the dynamics of vehicular traffic.

A network model called the ribosome flow model network with a finite pool (RFMNP) includes several RFMs that are interconnected to each other via a single pool of free particles [25]. It analyses the behaviour of such networks where the entry rate of the particles on lanes is affected by a single pool, that is, a single pool feeds all the input to the RFMs, and the output of each RFM is fed again into the pool. The number of particles remains conserved in this network and the network properties have been rigorously studied using tools from cooperative dynamical systems with a first integral [52–54]. An RFM network in [26] is a generalized network that analyses various network topologies using a set of interconnected RFMs. It models the static connections between the RFMs, and hence, the input to each RFM is a source (maybe pools of free ribosomes in the cell) whose output rate is a constant or proportion of the output of other RFMs. Here, the pool supplies a constant input source, and hence, it does not take into account the competition effects on the network’s behaviour owing to finite resources.

Almost all prior research has provided an understanding of biological activities constrained by a single pool. However, taking into account the concept of multiple pools helps one to comprehend the participation of particles in the vicinity of their targets. For instance, in the context of introducing synthetic circuit genes, a network model called orthogonal ribosome flow model (ORFM) was introduced, where the ribosome pool has been divided by the use of orthogonal ribosomes, and the introduced genes are only translated by mutated ribosomes [55]. The concept of orthogonal ribosomes was used to increase the protein output by decoupling circuit genes from the host pool of ribosomes. The concept of multiple pools also provides a useful framework to model vehicular flow between different cities where each pool represents each city.

In this article, we study the idea that the entry rate of the particles on lanes is affected by the occupancy of the nearby pool. Studying a minimal model of a two-pool network is a useful strategy for gaining insights into the behaviour of more complex systems with multiple pools. We present our theoretical investigation of a two-pool network and then illustrate how one can generalize it to study a network with multiple pools. We introduce a new model called the ribosome flow model network with two pools (RFMNTP) that include several RFMs interconnected via two dynamical pools of free particles. It captures the feature that the particles located far away from RFMs will not impact the initiation rates of these RFMs. Therefore, each RFM can be associated with two pools of particles, one pool containing particles impacting the initiation rate of the RFM and the other pool receiving its output. The whole system being closed conserves the number of particles. This two-pool network then models vehicular flow between the two cities where each pool represents each city. Similarly, pedestrian flow involving the movement of people between two places is reminiscent of the fact that it can be studied by incorporating two pools in the network.

By using the theory of cooperative dynamical systems with a first integral, we prove that the RFMNTP admits a continuum of steady-state points [52–54]. Therefore, the same steady-state point is attained by any two solutions starting from initial conditions corresponding to an equal total number of particles in the network, and hence, the network can be analysed by the steady-state density profile. This theoretical analysis can also be easily extended to prove the stability results for multiple pools. Next, we study how a change in one RFM affects the dynamics of the network. The change in the steady state with respect to changing the rate of a site in a specific RFM, say R , can be any one of the following outcomes: (i) steady-state pool densities can both increase (decrease) simultaneously, (ii) a decrease in the steady-state pool density that is feeding R and an increase in the steady-state density of the other pool. The results hold for any set of parameters in the network.

The structure of this article is organized as follows. The next section reviews the RFM and introduces the network of RFMs with multiple pools. A two-pool network and our main mathematical results are described in §3. The idea of understanding the dynamics of a network with multiple pools is illustrated in §4. The final section summarizes and suggests some directions for further research. The proofs of the results are provided in the Appendix for ease of reading. The source code which simulates differential equations described in this article can be accessed on GitHub at the following link: https://github.com/aditi081/matlab-codes-Networks-of-RFMs-with-multiple-pools.git.

2. The mathematical framework

Since the building blocks of our network model are RFMs, we first review RFM briefly.

2.1. The ribosome flow model

The RFM is a continuous-time model for the unidirectional flow of particles along a lattice of n consecutive sites. The state variable xj(t) , where j∈{1,2,…,n} represents the normalized particle density at site j at time t , where xj(t)=0 [xj(t)=1] represents that site j is totally empty [full]. The RFM is characterized by n+1 positive rates λj , where j∈{0,1,…,n} . In particular, the entry [exit] rate of particles is controlled by λ0 [ λn ]. The rate of flow from site j to site j+1 is given by λjxj(1−xj+1) , that is, the ribosomes are more likely to flow from site j to site j+1 when xj is full and xj+1 is empty. The output rate of exit of particles is R(t):=λnxn(t) . Figure 1 shows the topology of the RFM. The RFM is given by the following system of n first-order ordinary differential equations (ODEs):

(2.1) x˙1=λ0(1−x1)−λ1x1(1−x2),x˙2=λ1x1(1−x2)−λ2x2(1−x3),⋮x˙n=λn−1xn−1(1−xn)−λnxn.

Figure 1. The RFM models unidirectional flow along a lattice of n sites. The parameter λi>0 represents the transition rate from site i to site i+1 , where λ0 ( λn ) represents the initiation (exit) rate. The state variable xi(t) represents the density of site i at time t . The output rate at any time t is R(t) .

The trajectories of the RFM evolve on the n -dimensional unit cube denoted C:=[0,1]n , and it admits a unique steady-state density e in Int( C ), where Int( C ) denotes the interior of C [29]. This mathematical analysis is based on the theory of cooperative dynamical systems [33,34]. In [40], it has been shown that there exists a spectral representation for the mapping from λ0 , λ1 , … , λn to steady-state e given by (n+2)×(n+2) Jacobi matrix

A:=[0λ0−1/20…0λ0−1/20λ1−1/2…00λ1−1/20…0⋮0…00λn−1/20…0λn−1/20].

Note that A is a non-negative, symmetric and irreducible matrix, and hence, it admits a unique maximal eigenvalue σ>0 , and the entries of a corresponding eigenvector ζ∈Rn+2 are all positive for all i∈{1,2,…,n+2} [56]. It has also been proved in [40] that the steady-state values of the RFM satisfy

(2.2) ej=ζj+2λj1/2σζj+1,  for  j=1,2,…,n,

and the steady-state output rate satisfies

(2.3) R=1σ2.

To build a network of interconnected RFMs, the first step is to extend the RFM into a single-input–single-output (SISO) system. This is done by adding a time-varying measurable and bounded function v:R+⟶R+ to the RFM [39], where R+=R+∪{0} . In the context of translation, the function v represents the flow of ribosomes into the mRNA from the cell environment. The equations describing the RFM with the input and the output are as follows:

(2.4) x˙1=λ0v(1−x1)−λ1x1(1−x2),x˙2=λ1x1(1−x2)−λ2x2(1−x3),⋮x˙n=λn−1xn−1(1−xn)−λnxn,

and R(t)=λnxn(t) .

The system given by equation (2.4) is a monotone control system [57]. It can also be seen that C is an invariant set of the dynamics, that is, any trajectory emanating from any point in C shall remain in it for all time t≥0 .

2.2. Network of ribosome flow models with multiple pools

Now, we provide some insights into the dynamics of the flow of particles on several lanes connected via multiple pools. Consider a network consisting of M pools and N RFMs. Each pool is connected to at least one RFM that receives its input from the pool and at least one other RFM that feeds its output to the pool. The particles that are not attached to any RFM are present in the pools. The topology of this network can be represented using a directed multigraph where each node represents the pool and the directed edges represent the RFMs that indicate the flow of particles from one pool to the other pool (see figure 2).

Figure 2. The graph representation of the network with multiple pools where each node (circle) represents the pool and the directed edges (dashed lines having arrows) represent a chain of sites on which particles undergo RFM dynamics. The arrow of the edge pointing to the pool represents that the output of the RFM is feeding the pool.

Let RFM #i be described by the tuple Ai:={ni,Gi,xji,λki} for i=1,2,…,N , j=1,2…,ni , k=0,1,…,ni , where ni is the dimension of RFM #i;Gi:R+⟶R+ is the input function of RFM #i;xji s are the state variables and λki s are the positive transition rates along RFM #i . Let Pool #j density be described by zj(t) for j=1,2,…,M , where zj represents the average number of particles in the pool. Assume that Pool #j is feeding the input to RFM #k , k∈I , where I is a subset of the set {1,2,…,N} and let RFM #k′ , k′∈I′ , where I′ is a subset of the set {1,2,…,N}∖I , is feeding its output to the Pool #j . Thus, the dynamics of RFM #k is described by the following ODEs:

(2.5) x˙1k=λ0kGk(zj)(1−x1k)−λ1kx1k(1−x2k),x˙2k=λ1kx1k(1−x2k)−λ2kx2k(1−x3k),⋮x˙nkk=λnk−1kxnk−1k(1−xnkk)−λnkkxnkk,

and the dynamics is described by the following balance equation for each Pool #j :

(2.6) z˙j=∑k′∈I′λnk′k′xnk′k′−∑k∈Iλ0kGk(zj)(1−x1k).

The entry rate of particles into the RFMs is modulated by the occupancy level of the nearby pool. Note that there is no direct link between the RFMs in the network, and the interconnections are via the pool of particles. The pool outflow functions Gk describe the likelihood that the particles will attach to the RFMs. In other words, these functions model the competition for particles between the RFMs. Therefore, RFMs having a more effective initiation rate λ0kGk(zj) have more influx of particles into them.

Each state variable xji represents the normalized particle density and Gi gives non-negative output. Therefore, the state space of the network is

(2.7) Ω=[0,1]n1×[0,1]n2…[0,1]nN×[0,∞)M.

Let

(2.8) F(t):=∑j=1Mzj(t)+∑i=1N∑j=1nixji(t)

describe the total occupancy of particles in the network at any time t . Since it is a closed system, F(t) is the first integral of the dynamics. Note that all the particles can be accommodated in any of the pools.

Analysing such networks requires information about interconnected RFMs and the pools, and therefore, a minimal model of a two-pool network provides a useful starting point for studying the behaviour of particle flow in a network with multiple pools. We now shall begin our study with a new model RFMNTP that considers two dynamic pools in the network. The proposed model considers the dynamics of particles on various tracks, wherein on some tracks, particles are recruited from one pool and return to the other pool and vice versa. This is a primary study of a system having two pools in the framework of a network of RFMs. The findings can be generalized to complex systems with multiple pools.

3. The ribosome flow model network with two pools

We consider a network model consisting of several RFMs and two finite pools: Pool I and Pool II, and this models the dynamics of the flow of particles on several lanes interconnected via two pools (see figure 3a ). We represent those RFMs, say m≥1 in number, whose input is received through Pool I, and output is supplied to Pool II as RFMXs. For the reverse case, the RFMs, say n≥1 in number, whose input is received through Pool II, and output is supplied to Pool I are referred to as RFMYs. We call this network RFMNTP (see figure 3b ). Let RFMX # i is described by the tuple Ei:={ℓi,Gi,xji,λki} for i=1,2,…,m , j=1,2…,ℓi , k=0,1,…,ℓi , where ℓi is the dimension of RFMX # i;Gi:R+⟶R+ is the i th input function; xji s are the state variables, and λki s are the positive transition rates along RFMX # i . Similarly, consider RFMY # i represented by the tuple Fi:={pi,Hi,yji,ηki} for i=1,2,…,n , j=1,2…,pi , k=0,1,…,pi , where pi is the dimension of RFMY # i;Hi:R+⟶R+ is the i th input function, yji s are the state variables and ηki s are the positive transition rates along RFMY # i . The Pool I and Pool II densities at time t are modelled by z1(t) and z2(t) , respectively.

Figure 3. ( a ) Topology of the network with two pools and five lanes: particles from Pool I (Pool II) traverse lanes 1,2 (3,4,5) and join the Pool II (Pool I) and then traverse lanes 3,4,5 (1,2) and again join Pool I (Pool II). Hence, Pool I (Pool II) supplies its input to lanes 1,2 (3,4,5) and receives its output from lanes 3,4,5 (1,2). In the case of vehicular traffic, particle/lane/pool represents the car/road/city. ( b )Topology of the RFMNTP: the m RFMXs receive their input from Pool I and supply their output to Pool II and the n RFMYs receive their input from Pool II and supply their output to Pool I.

Thus, RFMX # i dynamics is described by the following equations:

(3.1) x˙1i=λ0iGi(z1)(1−x1i)−λ1ix1i(1−x2i),x˙2i=λ1ix1i(1−x2i)−λ2ix2i(1−x3i),⋮x˙ℓii=λℓi−1ixℓi−1i(1−xℓii)−λℓiixℓii,

and its output rate of exit of particles is given by λℓiixℓii .

The dynamics of RFMY # i is described by the following equations:

(3.2) y˙1i=η0iHi(z2)(1−y1i)−η1iy1i(1−y2i),y˙2i=η1iy1i(1−y2i)−η2iy2i(1−y3i),⋮y˙pii=ηpi−1iypi−1i(1−ypii)−ηpiiypii,

and its output rate of exit of particles is given by ηpiiypii .

Pool I feeds all the RFMXs, and the output of each RFMY is supplied into Pool I, so the change in z1 is given by the following balance equation:

(3.3) z˙1=∑i=1nηpiiypii−∑i=1mλ0i Gi(z1)(1−x1i).

In addition, Pool II feeds all the RFMYs, and the output of each RFMX is supplied into Pool II, so the change in z2 is given by the following balance equation:

(3.4) z˙2=∑i=1mλℓiixℓii−∑i=1nη0i Hi(z2)(1−y1i).

It can be observed that if the pools are empty, then no particles can attach to the respective lanes, and as the pools become fuller, more particles can attach to the lanes. Therefore, these properties are satisfied by imposing the following assumptions on Gi and Hi : (i) Gi(0)=0 and Hi(0)=0 and (ii) Gi(z1) and Hi(z2) are continuous and strictly increasing functions of z1 and z2 , respectively. Let

(3.5) Q(t):=z1(t)+z2(t)∑i=1m∑j=1ℓixji(t)+∑i=1n∑j=1piyji(t)

describe the total occupancy of particles in the network at any time t. Since the RFMNTP is a closed system, Q(t) is the first integral of the dynamics, that is, Q(t)=Q(0) for all t≥0 . Note that both pool densities are bounded by Q(0) . Summing up, the RFMNTP is a dynamical system with s=∑i=1mℓi+∑i=1npi+2 state variables whose dynamics are given by equations (3.1)–(3.4). The next section rigorously analyses the mathematical aspect of the RFMNTP.

3.1. Dynamical properties of the ribosome flow model network with two pools

In this section, we shall describe the various dynamical properties related to the RFMNTP. Given two vectors u , v∈Rn , we define order relation u≪v if ui<vi for all i . Recall that every xji and yji represent normalized particle density, and the assumptions on Gi and Hi imply that the pool densities are always non-negative. Therefore, the state space of the RFMNTP is

(3.6) B=[0,1]ℓ1×[0,1]ℓ2×…[0,1]ℓm×[0,1]p1×[0,1]p2×…[0,1]pn×[0,∞)×[0,∞).

Let [x(t,a) y(t,a) z1(t,a) z2(t,a)]′ denote the solution of the RFMNTP at time t for the initial condition a∈B , where x and y are the vector consisting of all the state variables of RFMXs and RFMYs, respectively. For r≥0 , let Lr:={a∈B:∑i=1sai=r} , that is, Lr represents all states in B corresponding to the total occupancy of particles equal to r in the network. In other words, Lr is the r level set of first integral Q .

3.2. Invariance

The following result states that for any initial condition a∈B , the trajectory of the RFMNTP stays in B for all t≥0 and follows by analysing the equations of the RFMNTP.

Proposition 3.1 The state space B is an invariant set of the RFMNTP, that is, 0≤xji(t,a)≤1 , 0≤yji(t,a)≤1 , zi(t,a)∈[0,∞) for any t≥0 and any initial condition a∈B .

We shall now show that the proposed nonlinear system of differential equations is a cooperative system (all the non-diagonal entries of the Jacobian matrix are non-negative in the invariant state space B ). This property guarantees the monotonicity (order-preserving property) of the flow with respect to the partial ordering in the phase space [34]. Now, the Jacobian matrix J of the vector field of the RFMNTP is

J(x,y,z1,z2)=[X100000…0U100X20000…0U20⋮⋱⋮⋮⋮⋮⋮⋮00…Xm00…0Um000…0Y10…00V100…00Y2…00V2⋮⋮⋮⋮⋱⋮⋮⋮00…000…Yn0VnA1A2…AmB1B2…BnZ10C1C2…CmD1D2…Dn0Z2],

where Xi represents the Jacobian matrix of RFMX #i and is given by

[−λ0iGi(z1)−λ1i(1−x2i)λ1ix1i000λ1i(1−x2i)−λ1ix1i−λ2i(1−x3i)λ2ix2i00⋱000−λℓi−2ixℓi−2i−λℓi−1i(1−xℓii)λℓi−1ixℓi−1i000λℓi−1i(1−xℓii)−λℓi−1ixℓi−1i−λℓii],

and Yi represents the Jacobian matrix of RFMY #i and is given by

[−η0iHi(z2)−η1i(1−y2i)η1iy1i000η1i(1−y2i)−η1iy1i−η2i(1−y3i)η2iy2i00⋱000−ηpi−2iypi−2i−ηpi−1i(1−ypii)ηpi−1iypi−1i000ηpi−1i(1−ypii)−ηpi−1iypi−1i−ηpii],

Ai=[λ0iGi(z1)0…00] ,    Bi=[0…0 0 ηpii] ,   Ci=[0…0 0 λℓii] ,   Di=[η0iHi(z2) 0…0  0] , Ui=[λ0iGi′(z1)(1−x1i)0...0]′ , Vi=[η0iHi′(z2)(1−y1i)0…0]′ , Z1=−∑i=1mλ0iGi′(z1)(1−x1i) and Z2=−∑i=1nη0iHi′(z2)(1−y1i) . Clearly, by proposition 3.1, we get that the Jacobian matrix J is Metzler for any initial condition in B , and thus, the RFMNTP is a cooperative dynamical system.

3.3. Persistence

The next result proves that the property called persistence holds, which implies that any trajectory becomes uniformly separated from the boundary of B .

Proposition 3.2 For any δ>0 there exists ϵ>0 depending on δ with ϵ→0 as δ→0 such that ϵ≤xji(t,a)≤1−ϵ , for i∈{1,2,…,m} , j∈{1,2,…,ℓi} , ϵ≤yji(t,a)≤1−ϵ , for i∈{1,2,…,n} , j∈{1,2,…,pi} and ϵ≤zi(t,a) for i∈{1,2} for all a∈(B∖{0}) and all t≥δ .

In other words, all the densities are smaller than one and larger than zero, and the pool occupancies are strictly positive. This result guarantees that the Jacobian matrix of the dynamics becomes irreducible after an arbitrarily short time [56]. Thus, the RFMNTP is a cooperative irreducible system of ODEs. Next, we analyse the asymptotic behaviour of the RFMNTP.

3.4. Stability

The following result shows that each level set of the first integral has a unique intersection with the ordered set of fixed points.

Theorem 3.1 The RFMNTP admits a unique steady-state point in every level set Lr , say er , and for any initial condition a∈Lr , the trajectory converges to er . Furthermore, for any 0≤r1<r2 , we have er1≪er2 .

The above theorem implies that the rates λji and ηji , and the total density of particle r determine a unique steady-state point of the network. Moreover, the continuum of the steady-state points {er:r∈[0,∞)} is linearly ordered. Combining proposition 3.2 and theorem 3.1, it follows that for any r>0 , er∈Int(B) , that is, the steady-state profile will never include densities of RFMXs and RFMYs that are either zero or one and the Pool I, and the Pool II steady-state densities are always strictly positive. The following example demonstrates theorem 3.1.

Example 3.1 Consider an RFMNTP with m=1 RFMX and n=1 RFMY each with dimension 2 . Assume that λ01=0.8 , λ11=1 , λ21=1.2 , η01=1 , η11=2 , η21=1 , G1(z1)=tanh⁡(z1) and H1(z2)=z2 . By theorem 3.1, there exists a unique equilibrium point e in L2 , and after simulating the dynamical system, we have e=[0.35890.23020.19090.27630.60230.3414]′ . Figure 4a,b depicts trajectories of RFMNTP for initial conditions in the level set L2:[0.5 0.5 0.5 0.5 0 0]′ and [0 0 0 0 1 1]′ , respectively. It can be observed that each of these trajectories converges to e .

Figure 4. Trajectories of the RFMNTP in example 3.1: ( a ) for initial condition, [0.5 0.5 0.5 0.5 0 0]′ and ( b ) for initial condition, [0 0 0 0 1 1]′ .

For ease of notation, let e=[ex ey ez1 ez2]′ where ex:=[ex11ex21⋯exℓ11ex12ex22⋯exℓ22⋯ex1mex2m⋯exℓmm] and ey:=[ey11ey21⋯eyp11ey12ey22⋯eyp22⋯eyp1neyp2n⋯eypnn] denote the unique steady-state point of the RFMNTP in the level set Lr of Q . At steady-state, the change in all the state variables and the pool variables with respect to time t becomes zero and thereby equation (3.1) implies

(3.7) λ0iGi(ez1)(1−ex1i)=λ1iex1i(1−ex2i)=⋯=λℓiiexℓii.

Similarly, at steady state, equation (3.2) yields

(3.8) η0iHi(ez2)(1−ey1i)=η1iey1i(1−ey2i)=⋯=ηpiieypii.

Again at steady state, equation (3.3) implies that

(3.9) ∑i=1nηpiieypii=∑i=1mλ0i Gi(ez1)(1−ex1i).

We can also express the above equation as

(3.10) ∑i=1nηpiieypii=∑i=1mλℓiiexℓii.

Equation (3.10) implies that the total output of all the RFMXs is equal to the total output of all the RFMYs.

Next, we provide an analysis of how the spectral approach can obtain the steady-state e of the network without any numerical simulations of the dynamics. Consider an RFMNTP with m RFMXs where each RFMX has dimension ℓi , i=1,2,…,m and n RFMYs, where each RFMY has dimension pk , k=1,2,…,n . We also assume that the transition rates in RFMX # i are represented by λji and in RFMY # k by ηki . The input to RFMX # i is given by Gi and to RFMY # k by Hk . Consider the total density of particles in the network to be r .

The steady-state values of each RFMX # i satisfy

(3.11) exji=ζj+2iλjiσxiζj+1i,  j=1,2,…,ℓi,

where σxi is the Perron eigenvalue and ζi is the corresponding Perron eigenvector of the (ℓi+2)×(ℓi+2) Jacobi matrix Ai given by

Ai:=[0(λ0iGi(ez1))−1/20…0(λ0iGi(ez1))−1/20(λ1i)−1/2…00(λ1i)−1/20…0⋮0…00(λℓii)−1/20…0(λℓii)−1/20].

The steady-state values of each RFMY # k satisfy

(3.12) eyjk=ξj+2kηjkσykξj+1k,  j=1,2,…,pk,

where σyk is the Perron eigenvalue and ξk is the corresponding Perron eigenvector of the (pk+2)×(pk+2) Jacobi matrix Bk given by

Bk:=[0(η0kHk(ez2))−1/20…0(η0kHk(ez2))−1/20(η1k)−1/2…00(η1k)−1/20…0⋮0…00(ηpkk)−1/20…0(ηpkk)−1/20].

It follows from equation (3.5) that

(3.13) ez1+ez2+∑i=1m∑j=1ℓiexji(ez1)+∑k=1n∑j=1pkeyjk(ez2)=r.

Also, equation (3.10) implies

(3.14) ∑i=1nηpiieypii(ez2)=∑i=1mλℓiiexℓii(ez1).

Combining equations (3.13) and (3.14) gives the expression of ez1 and ez2 in terms of the number of particles r . Thus, the entire steady-state profile of the network can then be calculated by plugging the values of r in the expression. Note that this approach also allows one to obtain an expression of densities for any unknown transition parameter and we only need to plug the values to obtain the entire steady-state profile without any numerical simulations of the dynamics. However, one needs to choose stable algorithms to solve the system of nonlinear equations (3.13) and (3.14).

3.5. Entrainment

Many important dynamic processes are periodic, such as the cell-cycle division process, gene regulation, circadian rhythm, 24 h solar day and more [58–60]. Proper functioning often requires such processes to vary periodically within the same period. For example, a person’s lack of synchronization between day and night can have health consequences [61]. For nonlinear systems, a periodic input signal does not guarantee that the response of the system will also be periodic as their behaviour can be quasi-periodic or chaotic [62,63]. Therefore, a natural question is whether the RFMNTP synchronizes with the periodic excitations or not. To answer this question, we assume that some or all of the parameters in the RFMNTP are not constants but periodic and continuous functions of time with a common period T>0 and satisfy the condition 0<δ1<λji≤δ2 and 0<δ3<ηji≤δ4 . In this case, we call the network model the periodic RFMNTP (PRFMNTP). The next result shows that all the trajectories approach a periodic pattern with the same period T .

Theorem 3.2. Consider the PRFMNTP. Fix r>0 . Then a unique T -periodic function ϕr:R+⟶Int(B) exists and for any a∈Lr , the solution of the PRFMNTP converges to ϕr .

In particular, PRFMNTP entrains to the periodic excitations in the parameters. As an additional point, if we examine the RFMNTP model for vehicular traffic between two cities, it enables a continuous flow of traffic while coordinating with the traffic lights. In simpler terms, when the traffic lights (rates) change periodically, the traffic density (state variables) will gradually settle into a recurring pattern within the same period. The following example illustrates the dynamic behaviour of the PRFMNTP model.

Example 3.2. Consider a PRFMNTP with m=1 RFMX with dimension ℓ1=3 and n=1 RFMY with dimensions p1=2 . Assume that λ01=0.8 , λ11=3 , λ21=3+2sin⁡(2πt) , λ31=3−2sin⁡(2πt) , η01=1.2 , η11=4−2sin⁡(2πt) , η21=1 , G1(z1)=z1 and H1(z2)=z2 . Let initial condition be xji=0 , yji=0 , z1(0)=0 and z2(0)=1 . Note that all the parameters are periodic with a common period T=1 . It can be seen in figure 5 that every trajectory converges to a periodic function with period one.

Figure 5. Trajectories of the PRFMNTP in example 3.2 as a function of time. Each state variable converges to a periodic solution having period one.

3.6. Effect of parameters

In this subsection, we analyse the effect of change in parameters on the network using a theoretical framework that is well explained through simple examples. Without loss of generality, we analyse the effect of change in the transition rate at a single site of an RFMX on the steady-state point of the network and presume that the change is in one of the transition rates in RFMX #1 . Let λ:=[λ01⋯λℓ11 λ02⋯λℓ22⋯λ0m⋯λℓmm] and η:=[η01⋯ηp11η02⋯ηp22⋯η0n⋯ηpnn] .

Theorem 3.3. Consider an RFMNTP with m RFMXs having dimensions ℓ1,ℓ2,…,ℓn and n RFMYs having dimensions p1,p2,…,pn . Let P=[λ  η]′ denotes the set of all parameters of the RFMNTP. Fix r>0 and let e denotes the unique steady-state point of the RFMNTP in the level set Lr of Q . Pick k∈{0,1,…,ℓ1} and suppose that we modify λk1 to λ¯k1 with λk1<λ¯k1 . Let e¯ denotes the steady-state point in the new RFMNTP. Then

(3.15) e¯xk1<exk1  and  exj1<e¯xj1  for all  j∈{k+1,…,ℓ1}.

Also, either one of the following cases holds:

ez1<e¯z1 , ez2<e¯z2 , exji<e¯xji for all i∈{2,3,…,m} , j∈{1,2,…,ℓi} and eyji<e¯yji for all i∈{1,2,…,n} , j∈{1,2,…,pi} .

e¯z1<ez1 , ez2<e¯z2 , e¯xji<exji for all i∈{2,3,…,m} , j∈{1,2,…,ℓi} and eyji<e¯yji for all i∈{1,2,…,n} , j∈{1,2,…,pi} .

e¯z1<ez1 , e¯z2<ez2 , e¯xji<exji for all i∈{2,3,…,m} , j∈{1,2,…,ℓi} and e¯yji<eyji for all i∈{1,2,…,n} , j∈{1,2,…,pi} .

ez1=e¯z1 , ez2<e¯z2 , exji=e¯xji for all i∈{2,3,…,m} , j∈{1,2,…,ℓi} and eyji<e¯yji for all i∈{1,2,…,n} , j∈{1,2,…,pi} .

e¯z1<ez1 , ez2=e¯z2 , e¯xji<exji for all i∈{2,3,…,m} , j∈{1,2,…,ℓi} and eyji=e¯yji for all i∈{1,2,…,n} , j∈{1,2,…,pi} .

Clearly, the above theorem lists all the possible cases of the effect of modifying the transition rate λk1 on the steady-state densities of the remaining RFMXs, all the RFMYs, and the pools. However, it does not give any information on the modified steady-state densities in sites {1,…,k−1} of RFMX #1 . Note that the theorem also exhibits how the output rates in all the RFMXs and RFMYs change. The next example demonstrates that the scenario when modifying a slow site increases the output rates of all lanes.

Example 3.3. Consider an RFMNTP with m=1 RFMX with dimension ℓ1=10 and n=2 RFMYs with dimensions pi=5 for i=1,2 . Assume that λj1=1 for i=1 , j=1,2,…,ℓi , ηji=1 for i=1,2 and j=1,2,…,pi , G1(z1)=z1 and Hi(z2)=z2 for i=1,2 . Let initial point be xji=0 , yji=0 , z1(0)=0.2 and z2(0)=0.2 . We simulate the system until steady state for a range of values of λ51 . It can be seen in figure 6a that we have ez1<e¯z1 and ez2<e¯z2 . Note that when λ51 is small, it is the only bottleneck rate in the RFMX and increasing it allows more particles to traverse RFMX more quickly. Hence, this increases the output flow rate of the RFMX, and subsequently, the Pool II density increases. This further increases the output rates of all the RFMYs and thus increases the Pool I density.

Figure 6. The steady-state pool densities for various values of transition rate λ51 in RFMX #1 of the RFMNTP considered in ( a ) example 3.3, ( b ) example 3.4 and ( c ) example 3.5.

It has been previously reported that owing to environmental change, stress conditions or pathological conditions, there could be particle stalling leading to an increase in a traffic jam in a track, resulting in a decrement in output rates of other tracks [25]. However, in our model, resource sharing is based on the concept that the entry rate is impacted by the neighbouring particles, and this can lead to both effects: a decrease in output rates from some of the tracks and an increase in output rates from others. The next example demonstrates this.

Example 3.4. Consider an RFMNTP with m=1 RFMX with dimension ℓ1=20 and n=1 RFMY with dimension p1=10 . Assume that λj1=1 for i=1 , j=1,2,…,ℓi except λ71=0.1 , ηji=1 for i=1 and j=1,2,…,pi , G1(z1)=z1 and H1(z2)=z2 . Consider an initial condition xji=0 , yji=0 , z1(0)=4 and z2(0)=4 . We simulate the system until steady state for a range of values of λ51 . It can be seen in figure 6b that we have e¯z1<ez1 and ez2<e¯z2 . This can be explained as follows. Note that when λ71 is the bottleneck rate in the RFMX, increasing λ51 only generates more traffic jams along RFMX. This depletes the Pool I density. However, in this case, the number of particles increases on the RFMX, and thus, the output flow rate of the RFMX increases, and subsequently, the Pool II density increases.

The following example exhibits the scenario in which increasing any of the transition rates in a specific lane yields an increase in the output rate of this lane, and the output rates in the other lanes all decrease.

Example 3.5. Consider an RFMNTP with m=2 RFMX with dimensions ℓ1=10 , ℓ2=5 and n=2 RFMYs with dimensions pi=5 for i=1,2 . Assume that λji=1 for i=1,2 , j=1,2,…,ℓi , ηji=1 for i=1,2 and j=1,2,…,pi , Gi(z1)=z1 and Hi(z2)=z2 for i=1,2 . Consider an initial condition xji=0 , yji=0 , z1(0)=4 and z2(0)=4 . We simulate the system until steady state for a range of values of λ51 . It can be seen in figure 6c that we have e¯z1<ez1 and e¯z2<ez2 . This can be understood by the following explanation. Increasing λ51 leads to the formation of traffic jams along RFMX # 1 owing to the bottleneck rate λ71 . This depletes Pool I and decreases the output rate of RFMX # 2 . So, there is a trade-off between the output values of both RFMXs, that is, whether the rate of increment of the output of RFMX # 1 is higher than the rate of decrement of RFMX # 2 . Depending upon the parameters of the RFMNTP, it can be seen in figure 6c that Pool II density decreased owing to the overall decrease in the total output rates from both RFMXs.

The next result provides specific information for the case when there is a single RFMX in the network.

Corollary 3.1. Consider an RFMNTP with m=1 RFMX and n RFMYs. Pick k∈{0,1,…,ℓ1} and suppose that λk1 is changed to λ¯k1 with λk1<λ¯k1 . Then,

(3.16) e¯xk1<exk1  and  exj1<e¯xj1  for all  j∈{k+1,…,ℓ1},

(3.17) ez2<e¯z2,  and  eyji<e¯yji  for all  i∈{1,2,…,n},j∈{1,2,…,pi}.

The following result implies that we can study steady-state properties of a network of m identical RFMXs and m identical RFMYs by a much simpler network consisting of only a single RFMX and a single RFMY.

Proposition 3.3. Consider the following two RFMNTPs:

(a) An RFMNTP with m identical RFMXs each having length ℓ , rates λ0,λ1…,λℓ and m identical RFMYs each having length p , rates η0,η1…,ηp . Let the output function G of Pool I and H of Pool II be homogeneous functions of degree 1 . Let r>0 be the total density of particles in the network and e denote its steady-state point.

(b) An RFMNTP with a single RFMX of length ℓ , rates (mλ0),λ1…,λℓ and a single RFMY of length p , rates (mη0),η1…,ηp . Let the output function G of Pool I and H of Pool II be homogeneous functions of degree one. Let r/m be the total density of particles in the network and e~=[e~x1e~x2⋯e~xℓ e~y1 e~y2⋯e~yp  e~z1 e~z2] denotes its steady-state point.

Then, we have

(3.18) exji=e~xj  for all  i=1,2,…,m,j=1,2,…ℓ,

(3.19) eyji=e~yj for all  i=1,2,…,m,j=1,2,…p,

and

(3.20) ez1=me~z1 and ez2=me~z2.

3.7. Mapping of the RFMNP to RFMNTP

In this section, we show that the RFMNP is a special case of our model RFMNTP. The RFMNP has been used for analysing the competition of ribosomes in the translation process. It assumes that the ribosomes that are located far away will also impact the initiation rates of the mRNAs and, therefore, include several RFMs interconnected via a single pool of free particles. All the RFMs feed the pool and the pool feeds the entry locations in all the RFMs. In this section, we show that the model RFMNP can be studied by the model RFMNTP, that is, we can construct the model RFMNP through RFMNTP as illustrated in the next paragraph.

Consider an RFMNP having m RFMs with dimensions ℓi , rates λji , state variables xji , a pool with density z , and the first integral having value (1/2)Q(0) . Construct an RFMNTP having m RFMXs and m RFMYs with the assumption ℓi=pi for all i=1,2,…,m , λji=ηji for all i=1,2,…,m , j=0,1,…,ℓi , Gi=Hi for each i=1,2…,m and having the first integral Q(0) . We shall show that both Pool I and Pool II steady-state density value is the same, that is, ez1=ez2 . Suppose on the contrary ez1<ez2 . Note that Gi is well defined and strictly increasing function and this implies Gi(ez1)<Gi(ez2) for all i. Since each RFMY is a copy of an RFMX, therefore, we have exℓii<eyℓii for all i , and this implies that ∑i=1mλℓii(eyℓii−exℓii)>0 , which is a contradiction to equation (3.10) in our case. Hence, ez1=ez2 , which implies that the steady-state densities of each RFMX # i are the same as of each RFMY # i , respectively. Therefore, our system becomes equivalent to the given one-pool network RFMNP.

3.8. Monte Carlo simulations

It has been shown that RFM and Monte Carlo simulations (TASEP with parallel update rule) provide highly correlated predictions for a large set of parameters [28]. In this section, we compare the steady states of the RFMNTP with the Monte Carlo simulations. This supports the modelling of the network of RFMs with two pools.

We validate our model by performing Monte Carlo simulations with a parallel update scheme. Each site is occupied with at most one particle, and the particle advances to the consecutive site if it is time to hop. The hopping times between consecutive sites of lanes receiving inputs from Pool I is exponentially distributed with parameters λ1i,λ2i,…,λℓii , that is, the next hopping time at site j of lane i is t+e(λji) , where t is the current time and e(λji) is randomly generated from the exponential distribution with mean parameter λji and the hopping time for the particle to hop to site 1 of lane i is calculated as t+e(λ0iGi(z1)) , where z1 is the number of particles present in Pool I at time t . Similarly, the hopping times between consecutive sites of lanes receiving inputs from Pool II are exponentially distributed with parameters η1i,η2i,…,ηpii the hopping time for the particle to hop to site 1 of lane i is calculated as t+e(λ0iHi(z2)) , where z2 is the number of particles present in Pool II at time t . A simulation begins with an empty chain with all the particles distributed arbitrarily in the pools and continues for 107 time steps. After removing the initial 104 steps from the calculation, the steady-state density of each site is calculated as the number of time steps it was occupied divided by the overall simulation runtime. In the example below, we show that simulations match the model RFMNTP.

Example 3.6 Consider an RFMNTP with m=2 RFMXs having dimensions ℓ1=10 , ℓ2=15 , n=1 RFMY having dimension p1=15 , λ0i=1 , λj1=1+θj , where θj is a random variable uniformly distributed in the interval (0,1) , λj2=2+θj , η01=1 , ηj1=5+θj , Gi(z1)=z1 , Hi(z2)=z2 and first integral having value 7 . It can be observed in figure 7 that the steady-state density profile of the RFMNTP and the Monte Carlo simulations match well with each other.

Figure 7. The steady-state density as a function of the site number for RFMXs and RFMYs in the RFMNTP in example 3.6. Solid lines and symbols denote numerically simulated RFMNTP and Monte Carlo simulations, respectively.

4. Analysing a network with multiple pools

Consider a network consisting of M pools and N RFMs having interconnections via the pools. One can extend the analysis done in §3.1 to show that the network admits a continuum of steady-state points. The following example demonstrates the dynamic behaviour of the network with three pools.

Example 4.1. Consider a network of three RFMs, each with dimension 2 , having three pools. Suppose RFM #1 receives its input from Pool #1 and supplies its output to Pool #2 , RFM #2 receives its input from Pool #2 and supplies its output to Pool #3 and RFM #3 receives its input from Pool #3 and supplies its output to Pool #1 . Assume that λ01=0.8 , λ11=1 , λ21=2 , λ02=1 , λ12=1.2 , λ22=0.1 , λ03=0.1 , λ13=0.5 , λ23=1 , G1(z1)=tanh⁡(z1) , G2(z2)=tanh⁡(z2) and G3(z3)=z3 . There exists a unique equilibrium point e in L4 and after calculation we have e=[0.095140.045410.82490.90820.19980.09080.126130.57451.1350]′ . Figure 8a,b depicts trajectories of RFMNTP for initial conditions in the level set L4: [0.5 0.5 0.5 0.5 0.5 0.5 0 0 1]′ and [0 0 0 0 0 0 1.5 1.5 1]′ , respectively. It can be observed that each of these trajectories converges to e .

Figure 8. Trajectories of the network with three pools in example 4.1: ( a ) for initial condition [0.5 0.5 0.5 0.5 0.5 0.5 0 0 1]′ and ( b ) for initial condition [0 0 0 0 0 0 1.5 1.5 1]′ .

Next, we understand the effect of modifying the transition rate of a site in an RFM on the entire network by analysing the steady-state densities. Equation (2.5) at steady-state e is given as

(3.21) λ0kGk(ezj)(1−ex1k)=λ1kex1k(1−ex2k)=⋯=λnkkexnkk.

Equation (2.6) at steady-state e is given as

(3.22) ∑k′∈I′λnk′k′exnk′k′=∑k∈Iλ0kGk(ezj)(1−ex1k).

Rewriting the above equation, we get

(3.23) ∑k′∈I′λnk′k′exnk′k′=∑k∈Iλnkkexnkk.

Without loss of generality, we assume that Pool #1 is feeding RFM #1 , and there is an increment in the rate λk1 of RFM #1 . Let e¯ represent the steady state of the modified network. Then by arguing similarly as in the proof of theorem 3.3, we get the information that e¯xk1<exk1 and exj1<e¯xj1for all j∈{k+1,…,n1} . The steady-state densities of the RFMs and the pools associated directly with Pool #1 follow the cases mentioned in theorem 3.3 depending on the various parameter values. The change in the steady-state densities of the other pools depends on the total input it is receiving from the RFMs and can be analysed through equation (3.23). Also, the case when ez1<e¯z1 and e¯zj<ezj for any j∈{1,2,…,M} is not possible as argued in theorem 3.3. This is a brief outline to gain an understanding of how the network with multiple pools behaves, as the exact scenario will be more clear when we know the interconnections between the RFMs via the pools.

In order to verify that the high correlation between the model and Monte Carlo simulations holds for a large set of parameters, we ran 250 tests, where in each test, a new set of rates are drawn randomly.

Example 4.2. Consider a network with M=3 pools and N=3 RFMs having dimensions n1=20 , n2=30 , n3=40 , where RFM #1 /RFM #2 /RFM #3 receives its input from Pool #1 /Pool #2 /Pool #3 and supplies its output to Pool #2 /Pool #3 /Pool #1 . Assume that λ0i=1 , λj1=0.5+θj , λj2=2+θj , λj3=1+θj , where θj is a random variable uniformly distributed in the interval (0,1) , Gi(zj)=zj and first integral having value 6 . Figure 9 depicts the correlations between the steady-state mean densities ( ρ ) of the RFMs and the steady-state mean densities ( σ ) calculated through Monte Carlo simulations. It can be seen that the correlation between the two is high ( r≃0.919633 ).

Figure 9. Steady-state mean densities (numerically simulated ρ and Monte Carlo simulated σ ) and the corresponding Pearson’s correlation coefficient r and p -value in example 4.2.

5. Discussion

Various transport phenomena involve the movement of particles along some tracks, for example, there is movement of RNA polymerases along DNA molecules, movement of ribosomes along mRNA molecules, motor proteins move along microtubules in order to transport cargo from one location to another, data packets move along buffers, there is the movement of vehicles along roads and many more. In all these systems, the particles are moving on a network of interconnected lanes in several transport systems. A common attribute in such phenomena is the presence of finite resources generating a closed network. The RFM was developed for analysing the excluded flow of particles along a one-dimensional isolated track, and it provides a useful and versatile modelling component that helps to understand the complex networks of the cellular as well as physical processes.

The network models consisting of a single pool are used to describe the behaviour of the system when the particles are distributed uniformly throughout the system. These single-pool models, however, do not take into account the distribution of particles in a local neighbourhood and, hence, are not able to model the movement of resources between different pools, for instance, the movement of cars (resources) between the two cities (pools). We introduced a new network model, RFMNTP, composed of several RFMs that focused on analysing how the network behaves when only the nearby resources impact the entry rates of its target. The RFMNTP is a closed network consisting of RFMs strategically connected to two pools such that Pool I (Pool II) feeds the input of some of the RFMs (remaining RFMs), and the output of them is fed into Pool II (Pool I). In other words, the first sites of some of the lanes and the last sites of the remaining lanes are connected to the same pool.

Understanding the stability of a system is a fundamental and foremost aspect of analysing systems in various fields of study. In this context, it is important to understand the stability of our network to predict its long-term behaviour. We prove that the RFMNTP is a cooperative irreducible dynamical system that admits a non-trivial first integral and, thus, enjoys several dynamical properties. In particular, the RFMNTP admits a continuum of steady-state points, and it entrains to periodic excitations in the parameters. Our theoretical analysis shows that an increase in the transition rate of a site in an RFM has a non-trivial effect on the output rates of the other RFMs. It can lead to any of the scenarios: the output rate of all the other RFMs increase or decrease; an increase in output rates of some of the RFMs and a decrease in output rates of the other RFMs. The specific outcome can be predicted by simulating the RFMNTP.

A noteworthy observation is that there could be a simultaneous increase in the output rates of some of the RFMs and a decrease in the output rates of the other RFMs. In the previous network model [25], we have seen that an increase in a transition rate in an RFM in the presence of a bottleneck rate leads to a decrease in the output rate of the other RFMs, whereas we can see in example 3.4 that this may not hold owing to aspect of local sharing of particles incorporated through two pools. Next, we have illustrated how to gain an understanding of how the changes in an RFM affect other RFMs and the overall behaviour of the network with multiple pools.

The model described here can be generalized to capture more complicated features. For example, the output of the shorter lanes can be fed back into the same pool. This phenomenon may be studied by adding its output rate to the same pool. The RFM also provides an analytical framework for modelling and analysing linear communication networks [64]. In this context, the moving particles are data packets, the chain of sites is a one-dimensional chain of ordered buffers, and the decreasing entry rate to a fuller buffer represents a kind of decentralized backpressure flow control. Another research direction is to analyse networks, comprising multiple flows that share common nodes, using a set of interconnected RFMs, constraining the link capacities in the communication networks. An applicability of our model can be to analyse a network topology where a common source node is linked to several chains of ordered buffers. The output of these chains at the destination node can be a source node for other sets of chains of ordered buffers and so on. One may also generalize RFMNTP by considering nearest-neighbour interactions in the network as seen in molecular motor traffic. Another interesting direction is to try and validate our predictions about the local behaviour of the cellular environment experimentally.

Ethics

This work did not require ethical approval from a human subject or animal welfare committee.

Data accessibility

This article has no additional data.

Declaration of AI use

We have not used AI-assisted technologies in creating this article.

Authors’ contributions

A.J.: conceptualization, formal analysis, investigation, methodology, software, validation, writing—original draft; A.K.G.: conceptualization, funding acquisition, methodology, project administration, supervision, writing—review and editing.

Both authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

We declare we have no competing interests.

Funding

No funding has been received for this article.

Appendix A. Proofs

Proof of Proposition 3.2: For simplicity, we consider RFMNTP with one RFMX and one RFMY. Let xj1:=xj and yj1:=yj . Also, let x0:=z1 , xℓ+1:=z2 , xℓ+j+1:=yj for j=1,2,…,p and x−1:=xℓ+p+1 . We now show that the system with state-variables x0,…,xℓ+p+1 satisfies the cyclic boundary-repelling (CBR) property given in [25] (i.e. for any δ>0 and any sufficiently small △>0 , ∃ P=P(δ,△)>0 such that for each k=1,…,ℓ+p+1 and each t≥0 the condition xk(t)≤△ and xk−1(t)≥δ implies that x˙k≥P ). For k=0 , we have x˙0=ηpxℓ+p+1−λ0G(x0)(1−x1)≥ηpδ−λ0G(△)(1−x1) , and we have G(0)=0 , G is a continuous function and thus x˙0≥ηpδ/2 for all △>0 sufficiently small. Now for k=1 , we have x˙1=λ0G(z1)(1−x1)−λ1x1(1−x2)≥λ0G(δ)(1−△)−λ1△(1−x2) , so x˙1≥λ0G(δ)/2 for all △>0 sufficiently small. For k∈{2,…,ℓ} , we have x˙k=λk−1xk−1(1−xk)−λkxk(1−xk+1)≥λkδ(1−△)−λk△(1−xk+1) , and therefore, x˙k≥λkδ/2 for all △>0 sufficiently small. For k=ℓ+1 , we have x˙ℓ+1=λℓxℓ−η0H(xℓ+1)(1−xℓ+2)≥λℓδ−η0H(△)(1−xℓ+2) , and we have H(0)=0 , H is a continuous function and thus x˙ℓ+1≥λℓδ/2 for all △>0 sufficiently small. Now, for k=ℓ+2 , we have x˙ℓ+2=η0H(z2)(1−xℓ+2)−η1xℓ+2(1−xℓ+3)≥η0H(δ)(1−△)−η1△(1−xℓ+3) , so x˙ℓ+2≥η0H(δ)/2 for all △>0 sufficiently small. For k∈{ℓ+3,…,ℓ+p+1} , we have x˙k=ηk−1xk−1(1−xk)−ηkxk(1−xk+1)≥ηkδ(1−△)−ηk△(1−xk+1) , and therefore, x˙k≥ηkδ/2 for all △>0 sufficiently small. Thus, the RFMNTP satisfies CBR.

Also, observe that if xk(t)>0 and t>0 then xk(T)>0 for all T≥t . It now follows from the result on repelling boundaries and persistence in ([25], lemma 1) that for any τ>0 , there exists ϵ1(τ)>0 , with ϵ1(τ)→0 as τ→0 , such that for any non-zero initial condition and any t≥τ the solution of the RFMNTP satisfies

(A 1) ϵ1≤xi(t), for all i∈{0,…,ℓ+p+1}.

For the other part of the equality, define ui(t):=1−xℓ+1−i(t) for i=1,2,…,ℓ and vi(t):=1−yp+1−i(t) for i=1,2,…,p . From equation (A 1), we have z1≥ϵ1 and z2≥ϵ1 for all t≥τ . Note that the system with state variables u1,u2,…,uℓ is an RFM having time-varying exit rate λ0G(z1) and the system with state variables v1,v2,…,vp is an RFM having time-varying exit rate η0H(z2) . Therefore, ∃ ϵ2(τ)>0 and ϵ3(τ)>0 such that any t≥τ , we have ϵ2≤ui(t), for all i∈{1,…,ℓ} and ϵ3≤vi(t), for all i∈{1,…,p} . Combining this with the definition of ui and vi completes the proof.∎

Proof of Theorem 3.1: We have that the RFMNTP is a cooperative irreducible system on Int(B) with a non-trivial first integral. Combining this with proposition 3.2 and the results in ([54], theorems 10 and 11) completes the proof of this theorem.∎

Proof of Theorem 3.2: From proposition 3.2, it follows that ϕr∈Int(B) , and the results in [53] (see theorem 3.1) prove this theorem.∎

Proof of Theorem 3.3: At steady-state, we have

(A 2) ∑i=1m∑j=1ℓiexji+∑i=1n∑j=1pieyji+ez1+ez2=r,

and this also holds for the modified network since the initial condition is the same, that is,

(A 3) ∑i=1m∑j=1ℓie¯xji+∑i=1n∑j=1pie¯yji+e¯z1+e¯z2=r.

Pick k∈{1,2,...ℓ1−1} . Let us assume that

(A 4) e¯xℓ11≤exℓ11.

Then equation (3.7) implies that

(A 5) λℓ1−11e¯xℓ1−11(1−e¯xℓ11)≤λℓ1−11exℓ1−11(1−exℓ11).

From equation (A 4), the above equation implies that e¯xℓ1−11≤exℓ1−11. Continuing this way, we have

(A 6) e¯xj1≤exj1  for all j=k+1,k+2,…,ℓj.

Now, from equation (3.7), consider

(A 7) λ¯k1e¯xk1(1−e¯xk+11)≤λk1exk1(1−exk+11).

We have λk1<λ¯k1 , thereby equation (A 7) implies that

(A 8) λ¯k1e¯xk1(1−e¯xk+11)<λ¯k1exk1(1−exk+11),

which implies e¯xk1<exk1 . From equation (3.7), we have

(A 9) λk−11e¯xk−11(1−e¯xk1)<λk−11exk−11(1−exk1).

Now, since e¯xk1<exk1 , we must have e¯xk−11<exk−11 . Continuing in this way, we get

(A 10) e¯xj1<exj1  for all j=1,2,…,k−2.

This also implies that e¯z1<ez1 . Since RFMX is a monotone control system, and therefore we have

(A 11) e¯xℓji<exℓji  for all i=2,3,…,m  and j=1,2,…,ℓi.

Note that all the RFMYs are interconnected through the pool variable z2 , and hence, all are affected in the same manner. From equations (3.10), (A 4) and (A 11), we have

(A 12) e¯yℓji<eyℓji  for all i=1,2,…,n  and j=1,2,…,pi,

and this implies e¯z2<ez2 , which yields the contradiction to the equation (A 3). Hence, exℓ11<e¯xℓ11 . Now if exk1≤e¯xk1 , this implies that ez1<e¯z1 which further implies eyℓji<e¯yℓji for all i=1,2,…,n and j=1,2,…,pi and ez2<e¯z2 . This is again a contradiction to equation (A 3). Also, note that the case ez1 is increasing and ez2 is decreasing is not possible as this will lead to the contradiction to equation (3.10). Now, consider k=0 . Let us assume that

(A 13) e¯xℓ11≤exℓ11.

Continuing as above, we have

(A 14) e¯xj1≤exj1  for all j=1,2,…,ℓ1−1.

Now, from equation (3.7), consider

(A 15) λ01G1(e¯z1)(1−e¯x11)≤λ01G1(ez1)(1−ex11).

We have λ01<λ¯01 and thereby equations (A 14) and (A 15) imply e¯z1<ez1 . Again using the above arguments, we get the contradiction to equation (A 13).

Now, consider k=ℓ1 . Then, we have to show that exℓ1−1<exℓ11 . Seeking a contradiction, assume

(A 16) exℓ11≤e¯xℓ11.

Then, equation (3.7) with the fact that λℓ11<λ¯ℓ11 implies that

(A 17) λℓ1−11exℓ1−11(1−exℓ11)<λℓ1−11e¯xℓ1−11(1−e¯xℓ11).

From equation (A 16), the above equation implies that exℓ1−11<e¯xℓ1−11. Continuing this way, we have

(A 18) exj1<e¯xj1  for all j=1,2,…,ℓ1−2.

This also implies that ez1<e¯z1 . Since RFMX is a monotone control system, and therefore, we have

(A 19) exℓji<e¯xℓji  for all i=2,3,…,m  and j=1,2,…,ℓi.

Note that all the RFMYs are interconnected through the pool variable z2 , and hence all are affected in the same manner. From equations (3.10), (A 16) and (A 19) and the fact that λℓ11<λ¯ℓ11 , we have

(A 20) eyℓji<e¯yℓji  for all i=1,2,…,n  and j=1,2,…,pi,

and this implies ez2<e¯z2 , which yields a contradiction to equation (A 3). Hence, e¯xℓ11<exℓ11 . This completes the proof of this theorem.∎

Proof of Proposition 3.3: The network considered in part (a) has identical RFMXs, and hence, the steady-state density profile of the RFMXs are the same, and similarly, this holds for RFMYs. Without loss of generality, consider the steady-state equations for RFMX #1 in (a)

(A 21) λ0G(ez1)(1−ex11)=λ1ex11(1−ex21)=⋯=λℓexℓ1,

and the steady-state equation for RFMY #1 in (a) is as follows:

(A 22) η0H(ez2)(1−ey11)=η1ey11(1−ey21)=⋯=ηpeyp1.

Similarly, the steady-state equation for RFMX in (b)

(A 23) mλ0G(e~z1)(1−e~x1)=λ1e~x1(1−e~x2)=⋯=λℓe~xℓ,

and the steady-state equation for the RFMY in (b) is as follows:

(A 24) mη0H(e~z2)(1−e~y1)=η1e~y1(1−e~y2)=⋯=ηpe~yp.

Now, consider equations (A 21) and (A 23); suppose we have

(A 25) λℓexℓ1=λℓe~xℓ

(A 26) ⟹exi1=e~xi for all i

and also

(A 27) λ0G(ez1)(1−ex11)=mλ0G(e~z1)(1−e~x1)

(A 28) ⟹ez1=me~z1.

Similarly, we get

(A 29) eyi1=e~yi  for all  i and ez2=me~z2.

Now, steady-state equation (3.5) for the network (a) is

(A 30) m(ex11+ex21+⋯+exℓ1)+m(ey11+ey21+⋯+eyp1)+ez1+ez2=r.

By equations (A 26), (A 28) and (A 29), we get

(A 31) e~x1+e~x2+⋯+e~xℓ+e~y1+e~y2+⋯+e~yp+e~z1+e~z2=r/m.

and hence, we get the steady-state equation for the network (b), and this completes the proof of the proposition.∎
==== Refs
References

1. Crick F . 1970 Central dogma of molecular biology. Nature 227 , 561–563. (10.1038/227561a0)4913914
2. Churchward G , Bremer H , Young R . 1982 Transcription in bacteria at different DNA concentrations. J. Bacteriol. 150 , 572–581. (10.1128/jb.150.2.572-581.1982)6175615
3. Vind J , Sørensen MA , Rasmussen MD , Pedersen S . 1993 Synthesis of proteins in Escherichia coli is limited by the concentration of free ribosomes: expression from reporter genes does not always reflect functional mRNA levels. J. Mol. Biol. 231 , 678–688. (10.1006/jmbi.1993.1319)7685825
4. Gyorgy A , Del Vecchio D . 2014 Limitations and trade-offs in gene expression due to competition for shared cellular resources. In IEEE 53rd Annual Conf. on Decision and Control (CDC), Los Angeles, CA, pp. 5431–5436. IEEE.
5. Frei T , Cella F , Tedeschi F , Gutiérrez J , Stan GB , Khammash M , Siciliano V . 2020 Characterization and mitigation of gene expression burden in mammalian cells. Nat. Commun. 11 , 4641. (10.1038/s41467-020-18392-x)32934213
6. Qian Y , Del Vecchio D . 2015 Effective interaction graphs arising from resource limitations in gene networks. In 2015 American Control Conf. (ACC), Chicago, IL, pp. 4417–4423 IEEE.https://ieeexplore.ieee.org/xpl/mostRecentIssue.jsp?punumber=7160954
7. Darlington APS , Kim J , Jiménez JI , Bates DG . 2018 Dynamic allocation of orthogonal ribosomes facilitates uncoupling of co-expressed genes. Nat. Commun. 9 , 695. (10.1038/s41467-018-02898-6)29449554
8. Ha M , Den Nijs M . 2002 Macroscopic car condensation in a parking garage. Phys. Rev. E. Stat. Nonlin. Soft Matter Phys. 66 , 036118. (10.1103/PhysRevE.66.036118)12366195
9. Wolf DE , Schreckenberg M , Bachem A . 1996 International workshop on traffic and granular flow. PO box 128, Farrer Road, Singapore: World scientific publishing co pvt ltd. See https://www.worldscientific.com/worldscibooks/10.1142/3077.
10. Chowdhury D , Santen L , Schadschneider A . 2000 Statistical physics of vehicular traffic and some related systems. Phys. Rep. 329 , 199–329. (10.1016/S0370-1573(99)00117-9)
11. Giridhar A , Kumar PR . 2006 Scheduling automated traffic on a network of roads. IEEE Trans. Veh. Technol. 55 , 1467–1474. (10.1109/TVT.2006.877472)
12. Jonathan Cook L , Zia RKP . 2012 Competition for finite resources. J. Stat. Mech. 2012 , 05008. (10.1088/1742-5468/2012/05/P05008)
13. Mather WH , Hasty J , Tsimring LS , Williams RJ . 2013 Translational cross talk in gene networks. Biophys. J. 104 , 2564–2572. (10.1016/j.bpj.2013.04.049)23746529
14. Sabi R , Tuller T . 2019 Modelling and measuring intracellular competition for finite resources during gene expression. J. R. Soc. Interface 16 , 20180887. (10.1098/rsif.2018.0887)31113334
15. Cook LJ , Zia RKP . 2009 Feedback and fluctuations in a totally asymmetric simple exclusion process with finite resources. J. Stat. Mech. 2009 , 2012. (10.1088/1742-5468/2009/02/P02012)
16. Greulich P , Ciandrini L , Allen RJ , Romano MC . 2012 Mixed population of competing totally asymmetric simple exclusion processes with a shared reservoir of particles. Phys. Rev. E 85 , 011142. (10.1103/PhysRevE.85.011142)
17. Neri I , Kern N , Parmeggiani A . 2011 Totally asymmetric simple exclusion process on networks. Phys. Rev. Lett. 107 , 068702. (10.1103/PhysRevLett.107.068702)21902376
18. MacDonald CT , Gibbs JH , Pipkin AC . 1968 Kinetics of biopolymerization on nucleic acid templates. Biopolymers 6 , 1–5. (10.1002/bip.1968.360060102)5641411
19. Zia RKP , Dong JJ , Schmittmann B . 2011 Modeling translation in protein synthesis with TASEP: a tutorial and recent developments. J. Stat. Phys. 144 , 405–428. (10.1007/s10955-011-0183-1)
20. Chowdhury D , Schadschneider A , Nishinari K . 2005 Physics of transport and traffic phenomena in biology: from molecular motors and cells to organisms. Phys. Life Rev. 2 , 318–352. (10.1016/j.plrev.2005.09.001)
21. Chou T , Mallick K , Zia RKP . 2011 Non-equilibrium statistical mechanics: from a paradigmatic model to biological transport. Rep. Prog. Phys. 74 , 116601. (10.1088/0034-4885/74/11/116601)
22. Kolomeisky AB . 1998 Asymmetric simple exclusion model with local inhomogeneity. J. Phys. A: Math. Gen. 31 , 1153–1164. (10.1088/0305-4470/31/4/006)
23. Chou T , Lakatos G . 2004 Clustered bottlenecks in mRNA translation and protein synthesis. Phys. Rev. Lett. 93 , 198101. (10.1103/PhysRevLett.93.198101)15600884
24. Dong JJ , Schmittmann B , Zia RKP . 2007 Towards a model for protein production rates. J. Stat. Phys. 128 , 21–34. (10.1007/s10955-006-9134-7)
25. Raveh A , Margaliot M , Sontag ED , Tuller T . 2016 A model for competition for ribosomes in the cell. J. R. Soc. Interface 13 , 116. (10.1098/rsif.2015.1062)
26. Nanikashvili I , Zarai Y , Ovseevich A , Tuller T , Margaliot M . 2019 Networks of ribosome flow models for modeling and analyzing intracellular traffic. Sci. Rep. 9 , 1703. (10.1038/s41598-018-37864-1)30737417
27. Jain A , Margaliot M , Gupta AK . 2022 Large-scale mRNA translation and the intricate effects of competition for the finite pool of ribosomes. J. R. Soc. Interface 19 , 20220033. (10.1098/rsif.2022.0033)35259953
28. Reuveni S , Meilijson I , Kupiec M , Ruppin E , Tuller T . 2011 Genome-scale analysis of translation elongation with a ribosome flow model. PLoS Comput. Biol. 7 , e1002127. (10.1371/journal.pcbi.1002127)21909250
29. Margaliot M , Tuller T . 2012 Stability analysis of the ribosome flow model. IEEE/ACM Trans. Comput. Biol. Bioinform. 9 , 1545–1552. (10.1109/TCBB.2012.88)22732691
30. Lohmiller W , Slotine JJE . 1998 On contraction analysis for non-linear systems. Automatica 34 , 683–696. (10.1016/S0005-1098(98)00019-3)
31. Russo G , di Bernardo M , Sontag ED . 2010 Global entrainment of transcriptional systems to periodic inputs. PLoS Comput. Biol. 6 , e1000739. (10.1371/journal.pcbi.1000739)20418962
32. Margaliot M , Sontag ED , Tuller T . 2016 Contraction after small transients. Automatica 67 , 178–184. (10.1016/j.automatica.2016.01.018)
33. Smillie J . 1984 Competitive and cooperative tridiagonal systems of differential equations. SIAM J. Math. Anal. 15 , 530–534. (10.1137/0515040)
34. Smith HL . 2008 Monotone dynamical systems: an introduction to the theory of competitive and cooperativesystems: an introduction to the theory of competitive and cooperative systems. Providence, Rhode Island, USA: AMS books. (10.1090/surv/041) See http://www.ams.org/surv/041
35. Margaliot M , Tuller T . 2012 On the steady-state distribution in the homogeneous ribosome flow model. IEEE/ACM Trans. Comput. Biol. Bioinform. 9 , 1724–1736. (10.1109/TCBB.2012.120)23221086
36. Kolomeisky AB , Schütz GM , Kolomeisky EB , Straley JP . 1998 Phase diagram of one-dimensional driven lattice gases with open boundaries. J. Phys. A: Math. Gen. 31 , 6911–6919. (10.1088/0305-4470/31/33/003)
37. Blythe RA , Evans MR . 2007 Nonequilibrium steady states of matrix-product form: a solver’s guide. J. Phys. A Math. Theor. 40 , R333–R441. (10.1088/1751-8113/40/46/R01)
38. Cook LJ , Zia RKP , Schmittmann B . 2009 Competition between multiple totally asymmetric simple exclusion processes for a finite pool of resources. Phys. Rev. E 80 , 031142. (10.1103/PhysRevE.80.031142)
39. Margaliot M , Tuller T . 2013 Ribosome flow model with positive feedback. J. R. Soc. Interface 10 , 20130267. (10.1098/rsif.2013.0267)23720534
40. Poker G , Zarai Y , Margaliot M , Tuller T . 2014 Maximizing protein translation rate in the non-homogeneous ribosome flow model: a convex optimization approach. J. R. Soc. Interface 11 , 20140713. (10.1098/rsif.2014.0713)25232050
41. Zarai Y , Margaliot M , Tuller T . 2017 A deterministic mathematical model for bidirectional excluded flow with Langmuir kinetics. PLoS ONE 12 , e0182178. (10.1371/journal.pone.0182178)28832591
42. Zarai Y , Margaliot M , Tuller T . 2017 Ribosome flow model with extended objects. J. R. Soc. Interface 14 , 135. (10.1098/rsif.2017.0128)
43. Margaliot M , Huleihel W , Tuller T . 2021 Variability in mRNA translation: a random matrix theory approach. Sci. Rep. 11 , 5300. (10.1038/s41598-021-84738-0)33674667
44. Jain A , Gupta AK . 2023 Modeling mRNA translation with ribosome abortions. IEEE/ACM Trans. Comput. Biol. Bioinform. 20 , 1600–1605. (10.1109/TCBB.2022.3203171)36044491
45. Edri S , Gazit E , Cohen E , Tuller T . 2014 The RNA polymerase flow model of gene transcription. IEEE Trans. Biomed. Circuits Syst. 8 , 54–64. (10.1109/TBCAS.2013.2290063)24681919
46. Zarai Y , Margaliot M , Kolomeisky AB . 2017 A deterministic model for one-dimensional excluded flow with local interactions. PLoS ONE 12 , e0182074. (10.1371/journal.pone.0182074)28796838
47. Jain A , Gupta AK . 2022 Modeling transport of extended interacting objects with drop-off phenomenon. PLoS ONE 17 , e0267858. (10.1371/journal.pone.0267858)35499998
48. Bar-Shalom E , Ovseevich A , Margaliot M . 2020 Ribosome flow model with different site sizes. SIAM J. Appl. Dyn. Syst. 19 , 541–576. (10.1137/19M1250571)
49. Herty M , Klar A , Singh AK . 2007 An ODE traffic network model. J. Comput. Appl. Math. 203 , 419–436. (10.1016/j.cam.2006.04.007)
50. Gupta AK , Redhu P . 2014 Analyses of the driver’s anticipation effect in a new lattice hydrodynamic traffic flow model with passing. Nonlinear Dyn. 76 , 1001–1011. (10.1007/s11071-013-1183-2)
51. Hossain MA , Tanimoto J . 2022 A microscopic traffic flow model for sharing information from a vehicle to vehicle by considering system time delay effect. Physica 585 , 126437. (10.1016/j.physa.2021.126437)
52. Mierczyński J . 1987 Strictly cooperative systems with a first integral. SIAM J. Math. Anal. 18 , 642–646. (10.1137/0518049)
53. Tang B , Kuang Y , Smith H . 1993 Strictly nonautonomous cooperative system with a first integral. SIAM J. Math. Anal. 24 , 1331–1339. (10.1137/0524076)
54. Mierczyński J . 2012 Cooperative irreducible systems of ordinary differential equations with first integral. arXiv 12084697.
55. Miller J , Al-Radhawi MA , Sontag ED . Mediating ribosomal competition by splitting pools. IEEE Control Syst. Lett. 5 , 1555–1560. (10.1109/LCSYS.2020.3041213)
56. Horn RA , Johnson CR . 2012 Matrix analysis. Cambridge, United Kingdom: Cambridge University Press. See https://www.cambridge.org/highereducation/product/9781139020411/book.
57. Angeli D , Sontag ED . 2003 Monotone control systems. IEEE Trans. Automat. Contr. 48 , 1684–1698. (10.1109/TAC.2003.817920)
58. Reppert SM , Weaver DR . 2002 Coordination of circadian timing in mammals. Nature 418 , 935–941. (10.1038/nature00965)12198538
59. Polymenis M , Aramayo R . 2015 Translate to divide: сontrol of the cell cycle by protein synthesis. Microb. Cell 2 , 94–104. (10.15698/mic2015.04.198)
60. Eser P , Demel C , Maier KC , Schwalb B , Pirkl N , Martin DE , Cramer P , Tresch A . 2014 Periodic mRNA synthesis and degradation co-operate during cell cycle gene expression. Mol. Syst. Biol. 10 , 717. (10.1002/msb.134886)24489117
61. Walker WH , Walton JC , DeVries AC , Nelson RJ . 2020 Circadian rhythm disruption and mental health. Transl. Psychiatry 10 , 28. (10.1038/s41398-020-0694-0)32066704
62. Kuznetsov YA . 2004 Numerical analysis of bifurcations. In Elements of applied bifurcation theory pp. 505–585. (10.1007/978-1-4757-3978-7)
63. Nikolaev EV , Rahi SJ , Sontag ED . 2018 Subharmonics and chaos in simple periodically forced biomolecular models. Biophys. J. 114 , 1232–1240. (10.1016/j.bpj.2018.01.006)29539408
64. Zarai Y , Mendel O , Margaliot M . 2015 Analyzing linear communication networks using the ribosome flow model. In 2015 IEEE Int. Conf. on Computer and Information Technology; Ubiquitous Computing and Communications; Dependable, Autonomic and Secure Computing; Pervasive Intelligence and Computing (CIT/IUCC/DASC/PICOM), Liverpool, UK, pp. 755–761 IEEE.
