The master equation plays an important role in many scientific fields including physics, chemistry, systems biology, physical finance and sociodynamics. We consider the master equation with periodic transition rates. This may represent an external periodic excitation like the 24 h solar day in biological systems or periodic traffic lights in a model of vehicular traffic. Using tools from systems and control theory, we prove that under mild technical conditions every solution of the master equation converges to a periodic solution with the same period as the rates. In other words, the master equation entrains (or phase locks) to periodic excitations. We describe two applications of our theoretical results to important models from statistical mechanics and epidemiology.

1. Introduction

Consider a physical system that can be in one of exactly N possible configurations and let x_i(t) denote the probability that the system is in configuration i at time t. We record the probabilities of all configurations at time t by the (column) state-vector $x (t) := [\begin{matrix} x_{1} (t) \\ ⋮ \\ x_{N} (t) \end{matrix}]$ . Every entry of this vector takes values in the closed interval [0,1].

The master equation describes the time evolution of these probabilities. It can be explained intuitively as describing the balance of probability currents going in and out of each possible state. To formulate the master equation for a specific model, one needs to know the rates of transition p_ij from configuration i to configuration j. A rigorous derivation of the master equation for a chemically reacting gas-phase system that is kept well stirred and in thermal equilibrium is given in [1]. The master equation plays a fundamental role in physics (where it is sometimes referred to as the Pauli master equation), chemistry, systems biology, sociodynamics and more. (For example, see the monographs [2,3] for more details.)

In this paper, we treat the general case where the transition rates p_ij may depend on both time t and on the probability distribution x(t) at time t. The resulting system of differential equations (see (1.2) below) constitutes a time-varying nonlinear dynamical system.

Note that even the special case, where the transition rates do not depend on the state x (the master equation (1.2) is then linear), is of general interest because it is intimately connected with the theory of Markov processes on finite configuration spaces. The relation is the following: such a Markov process is uniquely determined by an initial probability distribution on the configurations 1,…,N and by transition probabilities (Q_τ(t))_ij that denote the probabilities to be in configuration i at time t given that the system is in configuration j at time τ<t. These transition probabilities need to satisfy the Chapman–Kolmogorov equations which are essentially equivalent to the condition that the columns of the matrix Q_τ(t) satisfy the linear version of the master equation known as the forward equation in the theory of Markov processes [4,5].

In many physical systems, the number of possible configurations N can be very large. For example, the well-known totally asymmetric simple exclusion principle TASEP model (e.g. [6,7] and the references therein) includes a lattice of n consecutive sites, and each site can be either free or occupied by a particle, so the number of possible configurations is N=2ⁿ. In such cases, simulating the master equation and numerically calculating its steady state may be difficult even for small values of n and special methods must be applied (e.g. [7,8]).

Here, we are interested in deriving theoretical results that hold for any N. Specifically, we consider the case where the transition rates p_ij(t,x) are periodic in time t with a common period T>0. In this situation, we arrive at a T-periodic master equation. For such systems, we consider the problem of entrainment (or phase-locking):

Problem 1.1

Given a system described by a T-periodic master equation, determine if for every initial condition the probabilities x_i(t), i=1,…,N, converge to a periodic solution with period T. If this is so, determine if the periodic solution is unique or not.

In other words, if we view the transition rates as a T-periodic excitation, then the problem is to determine if the state of the system entrains, that is, converges to a periodic trajectory with the same period T. If this is so, an important question is whether there exists a unique periodic trajectory γ and then every solution converges to γ.

Entrainment is important in many natural and artificial systems. For example, organisms are often exposed to periodic excitations like the 24 h solar day and the periodic cell-cycle division process. Proper functioning often requires accurate entrainment of various biological processes to this excitation [9]. For example, cardiac arrhythmia is a heart disorder occurring when every other pulse generated by the sinoatrial node pacemaker is ineffective in driving the ventricular rhythm [10].

Epidemics of infectious diseases often correlate with seasonal changes and the required interventions, such as pulse vaccination, may also need to be periodic [11]. In mathematical population models, this means that the so-called transmission parameter is periodic, with a period of 1 year, and entrainment means that the spread of epidemics converges to a periodic pattern with the same period. As another example, traffic flow is often controlled by periodically varying traffic lights. In this context, entrainment means that the traffic flow converges to a periodic pattern with the same period as the traffic lights. This observation could be useful for the study of the green wave phenomenon [12]. Another example, from the field of power electronics, involves connecting a synchronous generator to the electric grid. The periodically varying voltage in the grid may be interpreted as a periodic excitation to the generator, and proper functioning requires the generator to entrain to this excitation (e.g. [13] and the references therein).

Our main results provide affirmative answers for problem 1.1 under quite general assumptions. Basic regularity assumptions on the transition probabilities that are required throughout the paper are summarized in assumption 2.1. In theorem 2.2, we then formulate with (2.5), see also (2.3), a condition that guarantees entrainment. It is observed in corollary 2.3 that condition (2.5) is always satisfied in the linear case. Uniqueness of the periodic attractor is shown in theorem 2.5 under the additional assumption of irreducibility, a condition that is well known in the theory of Markov processes (a definition of irreducibility is provided just before the statement of theorem 2.5).

In the special case of time-invariant rates, problem 1.1 reduces to determining if every solution converges to a steady state, and whether there exists a unique steady state. Indeed, time-invariant rates are T-periodic for any T>0 and thus entrainment means convergence to a solution that is T-periodic for any T>0, i.e. a steady state. This basic observation is the content of corollary 2.4.

As the x_is represent probabilities,

x_{i} (t) \in [0, 1] for all i and \sum_{j} x_{j} (t) = 1,

1.1

for all t≥t₀. The structure of the master equation guarantees that if x(t) satisfies (1.1) at time t=t₀, then (1.1) holds for all t≥t₀ even when the x_is are not necessarily linked to probabilities (see (1.9) below). Our results also hold of course in this case. The next example demonstrates such a case.

Example 1.2

An important topic in sociodynamics is the formation of large cities due to population migration. Haag [2], ch. 8 considers a master equation describing the flow of individuals between N settlements. The transition rates p_ij in this model represent the probability per time unit that an individual living in settlement i will migrate to settlement j. A mean-field approximation of this master equation yields a model in the form (1.2), where x_i represents the average density at settlement i, and $p_{i j} = \exp ((x_{j} - x_{i}) k_{i j})$ , with k_ij>0. This models the fact that the rate of transition from settlement i to settlement j increases when the population in settlement j is larger than in i, i.e. the tendency of individuals to migrate to larger cities. Note that the rates here are state-dependent, but not time-dependent. However, it is natural to assume that migration decisions depend on the season. For example, the tendency to migrate to colder cities may decrease (increase) in the winter (summer). This can be modelled by adding time dependence, say, changing the scaling parameters k_ij to functions k_ij(t) that are periodic with a period of 1 year. Then the transition rates depend on both state and time, and are periodic.

It is important to note that, in general, nonlinear dynamical systems do not entrain to periodic excitations. Indeed, Nikolaev et al. [14] discusses two ‘simple looking’ nonlinear dynamical systems whose response to periodic forcing is chaotic (rather than periodic). Moreover, these systems commonly appear as components of larger sensing and signal transduction pathways in systems biology. This highlights the importance of proving that entrainment does hold in specific classes of dynamical systems.

Although entrainment has attracted enormous research attention, it seems that it has not been addressed before for the general case of systems modelled using a T-periodic master equation. Here we apply the theory of cooperative dynamical systems admitting a first integral to derive conditions guaranteeing that the answer to problem 1.1 is affirmative. In §3, we describe two applications of our approach to important systems from statistical physics. The first is the totally asymmetric simple exclusion process (TASEP). This model has been introduced in the context of biocellular processes [15] and has become the standard model for the flow of ribosomes along the mRNA molecule during translation [16,17]. More generally, TASEP has become a paradigmatic model for the statistical mechanics of non-equilibrium systems [6,7,18]. It is in particular used to study the stochastic dynamics of interacting particle systems such as vehicular traffic [19].

The second application is to an important model from epidemiology called the stochastic susceptible–infected–susceptible (SIS) model.

The remainder of this paper is organized as follows. In the following subsection, we briefly explain the central mathematical concepts used in the proofs of our main results theorems 2.2 and 2.5. The exact mathematical formulation of these results is provided in §2. The subsequent section then describes the two applications to statistical physics and to epidemiology mentioned above. This is followed by a brief discussion of the significance of the results and an outlook on possible future directions of research in §4. The appendix includes all the proofs. These are based on known tools, yet we are able to use the special structure of the master equation to derive stronger results than those available in the literature on monotone dynamical systems.

1.1. Formulation of master equation and concepts of proof

We begin by formulating the master equation, determined by given transition rates p_ij(t,x)≥0, that governs the evolution of the probability distribution on N configurations:

\begin{aligned} {\dot{x}}_{1} (t) & = \sum_{\begin{matrix} j = 1 \\ j \neq 1 \end{matrix}}^{N} p_{j 1} (t, x (t)) x_{j} (t) - \sum_{\begin{matrix} j = 1 \\ j \neq 1 \end{matrix}}^{N} p_{1 j} (t, x (t)) x_{1} (t), \\ ⋮ \\ {\dot{x}}_{N} (t) & = \sum_{\begin{matrix} j = 1 \\ j \neq N \end{matrix}}^{N} p_{j N} (t, x (t)) x_{j} (t) - \sum_{\begin{matrix} j = 1 \\ j \neq N \end{matrix}}^{N} p_{N j} (t, x (t)) x_{N} (t) . \end{aligned}}

1.2

In lemma A.4 below, it is shown in a slightly more general setting that system (1.2) defines a flow in the set of probability distributions ${x \in {[0, 1]}^{N} ∣ \sum_{i = 1}^{N} x_{i} = 1}$ . For T>0, we refer to (1.2) as the T-periodic master equation if

p_{i j} (t + T, x) = p_{i j} (t, x),

1.3

for all i,j, all t and all x. Note that this includes the case where one (or several) of the rates is (are) T-periodic with T>0, and the other rates are time-independent, as a time-independent function satisfies (1.3) for all T. Clearly, from (1.3), it also follows that p_ij(t+kT,x)=p_ij(t,x) for any integer k. To make the period a well-defined notion one, therefore, often requires that the period is the minimal real number T>0 for which (1.3) is satisfied. Then constant functions do not have a period. As we want to include here the case of time-independent transition rates, e.g. corollary 2.3, we do not require the minimality of the common period T in (1.3).

For small values of N, it is sometimes possible to solve the master equation and then analyse entrainment directly. The next example demonstrates this.

Example 1.3

Consider the master equation (1.2) with N=2 and continuous time- (but not state-) dependent rates, i.e. p_ij=p_ij(t)≥0. Then (1.2) can be written as

[\begin{matrix} {\dot{x}}_{1} \\ {\dot{x}}_{2} \end{matrix}] = [\begin{matrix} - p_{12} & p_{21} \\ p_{12} & - p_{21} \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \end{matrix}] .

1.4

Assume also that all the rates are periodic with period T>0. Using the fact that x₁(t)+x₂(t)≡1 yields

{\dot{x}}_{1} (t) = p_{21} (t) - (p_{12} (t) + p_{21} (t)) x_{1} (t) .

1.5

Recall that x₁(0),x₂(0)∈[0,1] with x₁(0)+x₂(0)=1. Equation (1.5) implies that x₁(t)∈[0,1], for all t≥0, and thus x₂(t)∈[0,1], for all t≥0. Solving (1.5) yields

\begin{aligned} x_{1} (t) & = \exp (- \int_{0}^{t} (p_{12} (s) + p_{21} (s)) d s) (x_{1} (0) + c (t)) \\ and x_{2} (t) & = 1 - x_{1} (t), \end{aligned}}

1.6

where

c (t) := \int_{0}^{t} p_{21} (τ) \exp (\int_{0}^{τ} (p_{12} (s) + p_{21} (s)) d s) d τ

. To analyse if every solution converges to a periodic solution, we consider two cases.

Case 1: If p₁₂(t)+p₂₁(t)≡0, then (1.6) yields x(t)≡x(0), i.e. every point in the state-space is an equilibrium point. This means in particular that every solution is a periodic solution with period T.

Case 2: Assume that there exists a time t*∈[0,T) such that p₁₂(t*)+p₂₁(t*)>0. (Note that by continuity this in fact holds on a time interval that includes t*.) The solution (1.6) is periodic with period T if and only if x₁(T)=x₁(0), i.e. if and only if

x_{1} (0) = \frac{\exp (- \int_{0}^{T} (p_{12} (s) + p_{21} (s)) d s) c (T)}{1 - \exp (- \int_{0}^{T} (p_{12} (s) + p_{21} (s)) d s)} .

1.7

It is straightforward to show that the right-hand side in this equation is in [0,1], so in this case, there exists a unique periodic trajectory γ(t), with γ₁(0) equal to the expression in (1.7) and γ₂(0)=1−γ₁(0). To determine if every trajectory converges to γ, let z(t):=x(t)−γ(t), that is, the difference between the solution emanating from x₀ and the unique periodic solution. Then

{\dot{z}}_{1} = - (p_{12} + p_{21}) z_{1} and z_{1} (0) = x_{1} (0) - γ_{1} (0) .

As p₁₂(t)+p₂₁(t) is non-negative for all t, positive on a time interval and T-periodic, z₁(t) converges to zero and we conclude that any trajectory of the system converges to the unique periodic solution γ.

Of course, when N>2 and the rates depend on both t and x, this type of explicit analysis is impossible, and the proof of entrainment requires a different approach.

In general, proving that a time-varying nonlinear dynamical system entrains to periodic excitations is non-trivial. Rigorous proofs are known for two classes of dynamical systems: contractive systems and monotone systems with additional structure like a tridiagonal Jacobian [20] or admitting a first integral.

A system is called contractive if any two trajectories approach one another at an exponential rate [21,22]. Such systems entrain to periodic excitations [9,23]. An important special case is asymptotically stable linear systems with an additive periodic input u, that is, systems in the form

\dot{x} = A x + B u,

1.8

with

x \in R^{N}

A \in R^{N \times N}

a Hurwitz matrix,¹

u \in R^{M}

and

B \in R^{N \times M}

. In this case, x(t) converges to a periodic solution γ(t) and it is also possible to obtain a closed-form description of γ using the transfer function of the linear system [24]. We note that even in the case that the p_ijs in (1.2) do not depend on x, i.e. when (1.2) is linear in x, the master equation is not of the form (1.8) because the periodic influence in (1.2) enters through the transition rates p_ij and not through an additive input channel.

Next, we turn to the notion of a first integral. Define $H : R^{N} \to R$ by H(y):=y₁+⋯+y_N. Equation (1.2) implies that

\sum_{i = 1}^{N} {\dot{x}}_{i} (t) \equiv 0,

1.9

so that the value of H(x(t)) remains constant under the flow, that is, H is a first integral of (1.2).

A system is called monotone if its flow preserves a partial order, induced by an appropriate cone K, between its initial conditions [25]. An important special case of monotone systems is cooperative systems for which the cone K is the positive orthant. To explain this, define a partial ordering between vectors $a, b \in R^{n}$ by a≤b if every entry of a is smaller or equal to the corresponding entry of b. For example, for vectors in $R^{3}$

[1.1 3.12 11] \leq [1.2 4 11],

but

[1.1 3.12 11] ≰ [1.2 0 11] .

A system

\dot{x} = f (x)

is called cooperative if for any two initial conditions a,b with a≤b the solutions satisfy x(t,a)≤x(t,b) for any time t≥0. In other words, the dynamics preserves the ordering between the initial conditions.

Cooperative systems that admit a first integral entrain to periodic excitations. It is interesting to note that proofs of this property often follow from contraction arguments [26].

The master equation (1.2) is, in general, not contractive, although as we will show in theorem A.8 below it is on the ‘verge of contraction’ with respect to the ℓ₁ vector norm (see [27] for some related considerations). However, (1.2) admits a first integral and is often a cooperative system (see theorem A.9 below). In particular, when the rates do not depend on the state, i.e. p_ij=p_ij(t), then (1.2) is always cooperative.

2. Main results

We begin by specifying the exact conditions on (1.2) that are assumed throughout. For a set S, let int(S) denote the interior of S. For any time t, x(t) is an N-dimensional column vector that includes the probabilities of all N possible configurations. The relevant state-space is thus

Ω := {y \in R^{N} | y_{i} \geq 0 for all i, and \sum_{i = 1}^{N} y_{i} = 1} .

For an initial time t₀≥0 and an initial condition x(t₀), let x(t;t₀,x(t₀)) denote the solution of (1.2) at time t≥t₀. For our purposes, it will be convenient to assume that the vector field associated with system (1.2) is not only defined on the set Ω, but on all the closed positive cones

R_{+}^{N} := {x \in R^{N} | x_{j} \geq 0 for all 1 \leq j \leq N} .

Throughout this paper, we assume that the following condition holds.

Assumption 2.1

There exists T>0 such that the transition rates p_ij(t,x) are: continuous and non-negative on $[0, T] \times R_{+}^{N}$ ; continuously differentiable with respect to x on $[0, T) \times int (R_{+}^{N})$ and the derivative admits a continuous extension onto $[0, T) \times R_{+}^{N}$ ; and are jointly periodic with period T, that is,

p_{i j} (t + T, x) = p_{i j} (t, x),

2.1

for all i,j, all

t \in [0, \infty)

and all

x \in R_{+}^{N}

Let relint(Ω) denote the relative interior of Ω, that is,

relint (Ω) = {y \in R^{N} | y_{i} > 0 for all i, and \sum_{i = 1}^{N} y_{i} = 1} .

Note that if the rates are only defined on x∈Ω, with partial derivatives with respect to x_j on relint(Ω) with continuous extensions to Ω, then they can be extended to

R_{+}^{N}

so that the conditions in assumption 2.1 hold. For example, by defining them to be constant on rays through the origin and multiplied by a cut-off function χ(|x|₁), where χ is a smooth function with a compact support in

[0, \infty)

, satisfying χ(s)=1 for s=1, and where |x|₁ denotes the ℓ₁-norm of x.

We now determine the conditions guaranteeing that (1.2) is a cooperative dynamical system. Note that (1.2) can be written as

\dot{x} (t) = f (t, x) := A (t, x (t)) x (t),

2.2

where

A \in R^{N \times N}

is a matrix with entries

a_{i j} (t, x (t)) := {\begin{cases} p_{j i} (t, x (t)), & if i \neq j, \\ - \sum_{\begin{matrix} k = 1 \\ k \neq i \end{matrix}}^{N} p_{i k} (t, x (t)), & if i = j . \end{cases}

2.3

The Jacobian of the vector field f is the N×N matrix

J (t, x) := \frac{\partial f (t, x)}{\partial x} = A (t, x) + B (t, x),

2.4

where B is the matrix with entries

b_{i j} := \sum_{k = 1}^{N} x_{k} (\partial a_{i k} / \partial x_{j})

. Recall that a matrix

M \in R^{n \times n}

is called Metzler if every off-diagonal entry of M is non-negative. It follows that if

p_{j i} (t, x) + \sum_{k = 1}^{N} x_{k} (\partial a_{i k} (t, x) / \partial x_{j}) \geq 0

, for all i≠j, all t≥t₀ and all x∈Ω then J(t,x) is Metzler for all t≥t₀ and all x∈Ω.

We can now state our first result.

Theorem 2.2

Suppose that

p_{j i} (t, x) + \sum_{k = 1}^{N} \frac{\partial a_{i k} (t, x)}{\partial x_{j}} x_{k} \geq 0, for all i \neq j, t \geq t_{0}, x \in Ω .

2.5

Then, for any t₀≥0 and any x(t₀)∈Ω the solution x(t;t₀,x(t₀)) of (1.2) converges to a periodic solution with period T.

If the rates depend on time, but not on the state, i.e. p_ij=p_ij(t) for all i,j, then the condition in theorem 2.2 always holds, and this yields the following result.

Corollary 2.3

Ifp_ij=p_ij(t) for alli,j, then for anyt₀≥0 and anyx(t₀)∈Ωthe solutionx(t;t₀,x(t₀)) of (1.2) converges to a periodic solution with periodT.

Thus, theorem 2.2 describes a technical condition guaranteeing entrainment, and this condition automatically holds in the case where all the rates are functions of time only.

If the rates depend on the state, but not on time then we may apply theorem 2.2 for all T>0. Thus, the trajectories converge to a periodic solution with an arbitrary period, i.e. a steady state. This yields the following result.

Corollary 2.4

Ifp_ij=p_ij(x) for alli,jand in addition condition (2.5) holds, then for anyt₀≥0 and anyx(t₀)∈Ωthe solutionx(t;t₀,x(t₀)) of (1.2) converges to a steady state.

In some applications, it is useful to establish that all trajectories of (1.2) converge to a unique periodic trajectory. Recall that a matrix $M \in R^{n \times n}$ , with n≥2, is said to be reducible if there exists a permutation matrix P∈{0,1}^n×n, and an integer 1≤r≤n−1 such that $P^{'} M P = [\begin{matrix} B & C \\ 0 & D \end{matrix}]$ , where $B \in R^{r \times r}$ , $D \in R^{(n - r) \times (n - r)}$ , $C \in R^{r \times (n - r)}$ and $0 \in R^{(n - r) \times r}$ is a zero matrix. A matrix is called irreducible if it is not reducible. It is well known that a Metzler matrix M is irreducible if and only if the graph associated with the adjacency matrix of M is strongly connected [28], theorems 6.2.14 and 6.2.24.

Theorem 2.5

Suppose that the conditions in theorem2.2hold and, furthermore, that there exists a time t*≥t₀such that A(t*,x)+B(t*,x) is an irreducible matrix for all x∈Ω. Then (1.2) admits a unique periodic solution γ in Ω, with period T, and every solution x(t;t₀,x(t₀)) with x(t₀)∈Ω converges to γ at an exponential rate.

Example 2.6

Consider the system in example 1.3. This is of the form (2.2) with

A (t) = [\begin{matrix} - p_{12} (t) & p_{21} (t) \\ p_{12} (t) & - p_{21} (t) \end{matrix}] .

If there exists a time t* such that p₁₂(t*),p₂₁(t*)>0 then A(t*) is irreducible. We conclude that, in this case, all the conditions in theorem 2.5 hold, so the system admits a unique T-periodic solution γ and every trajectory converges to γ. This agrees of course with the results of the analysis in example 1.3 above where we arrived at the same conclusion under the slightly weaker assumption that p₁₂(t*)+p₂₁(t*)>0.

The next section describes an application of our results to two important models.

3. Applications

3.1. Entrainment in totally asymmetric simple exclusion process

The totally asymmetric simple exclusion process (TASEP) is a stochastic model of particles hopping along a one-dimensional chain. A particle at site k hops to site k+1 (the next site on the right) with an exponentially distributed probability² with rate h_k, provided the site k+1 is not occupied by another particle. This simple exclusion property generates an indirect link between the particles and allows to model the formation of traffic jams. Indeed, if a particle ‘gets stuck’ for a long time in the same site, then other particles accumulate behind it. At the left end of the chain particles enter with a certain entry rate α>0 and at the right end particles leave with a rate β>0 (figure 1).

Figure 1. The TASEP model includes particles randomly hopping along a chain of n sites. Note that the particle in site 1 cannot hop forward because site 2 contains a particle.

As pointed out in the introduction, TASEP has become a standard tool for modelling ribosome flow during translation, and is a paradigmatic model for the statistical mechanics of non-equilibrium systems. We note that in the classical TASEP model the rates α, β and h_i are constants, but several papers considered TASEP with periodic rates [29–31] that can be used, for example, as models for vehicular traffic controlled by periodically varying traffic signals.

It was shown in [32] that the dynamic mean-field approximation of TASEP, called the ribosome flow model (RFM), entrains. However, the RFM is not a master equation and the proof of entrainment in [32] is based on different ideas. For more on the analysis of the RFM, see e.g. [33–36].

For a chain of length n, denoting an occupied site by 1 and a free site by 0, the set of possible configurations is {0,1}ⁿ, and thus the number of possible configurations is N=2ⁿ. The dynamics of TASEP can be expressed as a master equation with transition rates p_ij that depend on the values α, β and h_i, i=1,…,n. For the sake of simplicity, we will show this in the specific case n=2, but all our results below hold for any value of n.

When n=2, the possible configurations of particles along the chain are C₁:=(0,0), C₂:=(0,1), C₃:=(1,0) and C₄:=(1,1). Let x_i(t) denote the probability that the system is in configuration C_i at time t, for example, x₁(t) is the probability that both sites are empty at time t. Then x₁ may decrease (increase) due to the transition C₁→C₃ [C₂→C₁], i.e. when a particle enters the first site (a particle in the second site hops out of the chain). This gives

{\dot{x}}_{1} (t) = - α x_{1} (t) + β x_{2} (t) .

Similar considerations for all configurations lead to the master equation $\dot{x} = A x$ , with

A := (\begin{matrix} - α & β & 0 & 0 \\ 0 & - α - β & h_{1} & 0 \\ α & 0 & - h_{1} & β \\ 0 & α & 0 & - β \end{matrix}) .

If the entry, exit and hopping rates are time dependent and periodic, all with the same period T, one easily sees that the resulting master equation satisfies assumption 2.1 as well as all assumptions of theorem 2.2. Hence, we conclude that every solution of the master equation starting in Ω converges to a periodic solution with period T. Moreover, if there exists a time t* such that α(t*),β(t*),h₁(t*)>0, then A(t*) is irreducible. Hence, the conditions of theorem 2.5 are also satisfied, so we conclude that the periodic solution is unique and convergence takes place at an exponential rate. It is not difficult to show that the same holds for TASEP with any length n.

Example 3.1

When n=3, the possible particle configurations are C₁:= (0,0,0), C₂:=(0,0,1), C₃:=(0,1,0), C₄:=(0,1,1),…,C₈:=(1,1,1). Let x_i(t) denote the probability that the system is in configuration C_i at time t. The TASEP master equation in this case is $\dot{x} = A x$ , with

A = [\begin{matrix} - α & β & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & - α - β & h_{2} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & - α - h_{2} & β & h_{1} & 0 & 0 & 0 \\ 0 & 0 & 0 & - α - β & 0 & h_{1} & 0 & 0 \\ α & 0 & 0 & 0 & - h_{1} & β & 0 & 0 \\ 0 & α & 0 & 0 & 0 & - h_{1} - β & h_{2} & 0 \\ 0 & 0 & α & 0 & 0 & 0 & - h_{2} & β \\ 0 & 0 & 0 & α & 0 & 0 & 0 & - β \end{matrix}] .

We simulated this system with the rates

α (t) = 1 + \cos (t), β (t) = 1 + \cos (t + π), h_{1} = \frac{1}{2}, h_{2} = \frac{1}{4}

and initial condition

x (0) = {[\frac{1}{8} \dots \frac{1}{8}]}^{'}

. Note that all the rates here are jointly periodic with period 2π. Figure 2 depicts x₁(t) (black square), x₄(t) (red asterisk) and x₈(t) (blue circle) as a function of t (we depict only three x_is to avoid cluttering the figure). Note that as the entry rate α(t) is maximal and the exit rate β(t) is minimal at t=0, the probability x₈(t) [x₁(t)] to be in state (1,1,1) [(0,0,0)] quickly increases (decreases) near t=0. As time progresses, the probabilities converge to a periodic pattern with period 2π.

Entrainment of the probabilities x_i has consequences for other quantities of interest in statistical mechanics. For instance, an important quantity is the occupation density, i.e. the probability that site k is occupied, often denoted by 〈τ_k〉, cf. [37,38]. Denoting the kth component of the configuration C_i∈{0,1}ⁿ by C_i,k, a straightforward computation reveals that

⟨ τ_{k} (t) ⟩ = \sum_{i = 1}^{N} C_{i, k} x_{i} (t) .

It is thus immediate that the occupation densities also converge to a unique periodic solution.

This phenomenon has already been observed empirically in [29] that studied a semi-infinite and finite TASEP coupled at the end to a reservoir with a periodic time-varying particle density. This models, for example, a traffic lane ending with a periodically varying traffic light. The simulations in [29] suggest that this leads to the development of a sawteeth density profile along the chain, and that ‘The sawteeth profile is changing with time, but it regains its shape after each complete period…’ [29], p. 011122-2 (see also [30,31] for some related considerations).

Our results can also be interpreted in terms of the particles along the chain in TASEP. As the expectation of the occupation densities 〈τ_k〉 converges to a periodic solution, this means that, in the long term, the TASEP dynamics ‘fluctuates’ around a periodic ‘mean’ solution (see e.g. the simulation results depicted in fig. 5 in [32]). Moreover, in [30,31] it was found for closely related models that the limiting periodic density profiles (whose existence is also guaranteed by our results) have an interesting structure that depends in a non-trivial way on the frequency of the transition rates.

3.2 Entrainment in a stochastic susceptible–infected–susceptible model

The stochastic SIS model plays an important role in mathematical epidemiology [39]. But, as noted in [40], it is usually studied under the assumption of fixed contact and recovery rates. Here, we apply our results to prove entrainment in an SIS model with periodic rates.

Consider a population of size N divided into susceptible and infected individuals. Let S(t) [I(t)] denote the size of the susceptible (infected) part of the population at time t, so that S(t)+I(t)≡N. We assume two mechanisms for infection. The first is by contact with an infected individual and depends on the contact rate a(t). The second is by some external agent (modelling, say, insect bite) with rate c(t). The recovery rate is b(t) (figure 3). We assume that a(t), b(t) and c(t) are continuous and take non-negative values for all time t.

Figure 3. Three mechanisms for transitions between the classes of susceptible and infected and the associated rates.

If I(t)=n (so S(t)=N−n), then the probability that one individual recovers in the time interval [t,t+dt] is b(t)n dt+o(dt), and the probability for one new infection to occur in this time interval is a(t)n((N−n)/N) dt+c(t)((N−n)/N) dt+o(dt). For n∈{0,…,N}, let P_n(t) denote the probability that I(t)=n. This yields the master equation:

{\dot{P}}_{n} = ((n - 1) a + c) (1 - \frac{n - 1}{N}) P_{n - 1} - ((n a + c) (1 - \frac{n}{N}) + n b) P_{n} + (n + 1) b P_{n + 1},

3.1

for n∈{0,1,…,N}, where we define P₋₁=P_N+1:=0, and for simplicity omit the dependence on t. This set of N+1 equations may be written in matrix form as

\dot{x} = M x,

where x:=[P₀P₁ … P_N]′∈[0,1]^N+1, M:=P−D, with D:=diag(0,b,2b,…,Nb), and P is the (N+1)×(N+1) matrix:

[\begin{matrix} - c q_{0} & b & 0 & 0 & \dots & 0 & 0 & 0 \\ c q_{0} & - (a + c) q_{1} & 2 b & 0 & \dots & 0 & 0 & 0 \\ 0 & (a + c) q_{1} & - (2 a + c) q_{2} & 3 b & \dots & 0 & 0 & 0 \\ 0 & 0 & (2 a + c) q_{2} & - (3 a + c) q_{3} & \dots & 0 & 0 & 0 \\ ⋮ \\ 0 & 0 & 0 & 0 & \dots & ((N - 2) a + c)) q_{N - 2} & - ((N - 1) a + c) q_{N - 1} & b N \\ 0 & 0 & 0 & 0 & \dots & 0 & ((N - 1) a + c) q_{N - 1} & 0 \end{matrix}],

where q_i:=1−i/N. Note that M(t) is Metzler, as a(t),b(t) and c(t) are non-negative for all t. Thus, theorems 2.2 and 2.5 yield the following result.

Corollary 3.2

Ifa(t), b(t) andc(t) are allT-periodic then any solution of (3.1) with $x (0) \in Ω \subset R^{N + 1}$ converges to aT-periodic solution. Furthermore, if there exists a timet*≥0 such that

b (t^{*}) c (t^{*}) > 0,

3.2

then there exists a uniqueT-periodic solutionγinΩand every solution converges toγ.

Example 3.3

Consider the stochastic SIS model with N=3, a(t)=1, $b (t) = 3 + 3 \cos (t + 0.5)$ and $c (t) = 2 - 2 \sin (t + 0.75)$ . These rates are non-negative and jointly T-periodic for T=2π and clearly there exists t*≥0 such that b(t*)c(t*)>0. Figure 4 depicts P_i(t), i=0,1,2, (note that P₃(t)=1−P₀(t)−P₁(t)−P₂(t)) as a function of time t for the initial condition $P (0) = (\frac{1}{4}) 1_{4} \in Ω$ . It may be seen that every P_i(t) converges to a periodic solution with period 2π. Taking other initial conditions x(0)∈Ω yields convergence to the same periodic solution.

Figure 4. Probabilities P₀(t) (black asterisk), P₁(t) (red circle) and P₂(t) (blue square) as a function of t in example 3.3.

Note that if the irreducibility condition (3.2) does not hold then the system may have several periodic solutions. To see this, consider, for example, the case b(t)=c(t)≡0. Let $e^{i} \in R^{N + 1}$ denote the vector with entry i equal to one and all other entries zero. Then both x(t)≡e¹ and x(t)≡e^N+1 are (periodic) solutions of the dynamics.

4. Discussion

In his 1929 paper on periodicity in disease prevalence, Soper [41] states: ‘Perhaps no events of human experience interest us so continuously, from generation to generation, as those which are, or seem to be, periodic.’ Soper also raised the question of whether the observed periodicity in epidemic outbreaks is the result of a ‘seasonal change in perturbing influences, such as might be brought about by school break-up and reassembling, or other annual recurrences?’ In modern terms, this amounts to asking whether the solutions of the system describing the dynamics of the epidemics entrain to periodic variations in the transmission parameters.

Here, we studied entrainment for dynamical systems described by a master equation. We considered a rather general formulation where the transition rates may depend on both time and state. Also, we did not assume any symmetry conditions (e.g. detailed balance conditions [3], ch. V) on the rates. We also note that this formulation implies similar results for nonlinear systems. Indeed, consider the time-varying nonlinear system:

\dot{x} = f (t, x)

4.1

and assume that f(t,0)=0, for all t. Let J(t,x):=(∂f/∂x)(t,x) denote the Jacobian of the vector field. Then

\begin{aligned} \dot{x} & = \int_{0}^{1} \frac{d}{d s} f (t, s x) d s \\ = A (t, x) x, \end{aligned}

where

A (t, x) := \int_{0}^{1} J (t, s x) d s

. If A(t,x) has the form (2.3) then the results above can be applied to (4.1).

We proved that entrainment indeed holds under quite mild technical conditions. This follows from the fact that the master equation is a cooperative dynamical system admitting a first integral. Owing to the prevalence of the master equation as a model for natural and artificial phenomena, we believe that this result will find many applications. To demonstrate this, we described two applications of our results: a proof of entrainment in TASEP and in a stochastic SIS model.

The rigorous proof that the solutions of the master equation entrain is of course a necessary first step in studying the structure of the periodic trajectory (or trajectories), and its dependence on various parameters. Indeed, in many applications it is of interest to obtain more information on the periodic trajectory, e.g. its amplitude. Of course, one cannot expect in general to obtain a closed-form description of the limit cycle. However, for contractive dynamical systems there do exist efficient methods for obtaining a closed-form approximation of the limit cycle accompanied by explicit error bounds [23]. Developing a similar approach for the attractive limit cycle of the master equation may be an interesting topic for further research. In the specific case of TASEP with fixed rates, there exists a powerful representation of the steady state in terms of a product of matrices [7,37]. It may be of interest to try and represent the periodic steady state using a similar product, but with matrices with periodic entries. This could be used in particular to study the effects of periodic perturbations to the boundary-induced phase transitions that have been observed for TASEP in [42].

Data accessibility

All the relevant data are contained in this paper.

Authors' contributions

M.M., L.G. and T.K. performed the research and wrote the paper together.

Competing interests

We declare we have no competing interests.

Funding

This research is supported, in part, by a research grant from the Israeli Science Foundation (ISF grant 410/15).

Acknowledgements

We thank Yoram Zarai for helpful comments. L.G. and T.K. are grateful to Joachim Krug for very helpful discussions on interacting particle systems and for pointing out a number of references.

Appendix A. Proofs of theorems 2.2 and 2.5

The proofs of theorems 2.2 and 2.5 are based on known tools from the theory of monotone dynamical systems admitting a first integral with a positive gradient (e.g. [43–45]). We present in this appendix a self-contained proof taking full advantage of the technical simplifications that our specific setting permits. This, in particular, allows us to prove that the results hold on the closed state-space and also that irreducibility at a single time point is enough to guarantee convergence to a unique periodic solution. Without loss of generality we always assume that the initial time is t₀=0. It is convenient to work with the ℓ₁ vector norm $| x |_{1} = \sum_{i} | x_{i} |$ .

We begin by introducing some notation. First recall the notation

R_{+}^{N} := {x \in R^{N} | x_{j} \geq 0 for all 1 \leq j \leq N},

for the closed positive cone. Define a set of vector fields by

F := {f : [0, \infty) \times R_{+}^{N} \to R^{N} | properties (i)–(v) below are satisfied},

where

(i) f is continuous;
(ii) for all j∈{1,…,N}, (∂f/∂x_j)(t,x) exists for $(t, x) \in [0, \infty) \times int (R_{+}^{N})$ and has a continuous extension onto $[0, \infty) \times R_{+}^{N}$ . Thus, J(t,x) from (2.4) is defined on $[0, \infty) \times R_{+}^{N}$ to be the continuous extension of (∂f/∂x_j)(t,x);
(iii) J(t,x) is Metzler for all $(t, x) \in [0, \infty) \times R_{+}^{N}$ ;
(iv) $\sum_{i = 1}^{N} f_{i} (t, x) = 0$ for all $(t, x) \in [0, \infty) \times R_{+}^{N}$ ;
(v) f(t,0)=0 for all $t \in [0, \infty)$ .

For T>0, let $F_{T} := {f \in F | f (t + T, x) = f (t, x) for all (t, x) \in [0, \infty) \times R_{+}^{N}}$ , that is, the set of vector fields in $F$ that are also T-periodic.

It is straightforward to check that f(t,x)=A(t,x)x with A defined by (2.3) belongs to $F_{T}$ if assumption 2.1 and the assumptions of theorem 2.2 hold. Therefore, theorem 2.2 follows from the following result.

Theorem A.1

If $f \in F_{T}$ then for all x₀∈Ω the solution x(t;x₀) of the initial value problem

\dot{x} = f (t, x), x (0) = x_{0},

is asymptotically T-periodic, i.e. there exists a solution

γ : R \to Ω

\dot{γ} = f (t, γ)

, with γ(t+T)=γ(t) for all

t \in R

, and

lim_{t \to \infty} | x (t; x_{0}) - γ (t) |_{1} = 0 .

Let

F_{irr}^{Ω} := {f \in F | there exists t^{*} \geq 0 such that J (t^{*}, x) is irreducible for all x \in Ω} .

The next result is a generalization of theorem 2.5.

Theorem A.2

If $f \in (F_{T} \cap F_{irr}^{Ω})$ then the differential equation $\dot{x} = f (t, x)$ admits a unique T-periodic solution $γ : R \to Ω$ . Moreover, there exists α>0 such that for any initial condition x₀∈Ω the corresponding solution x(t;x₀) satisfies

| x (t; x_{0}) - γ (t) |_{1} < 2 e^{- α t} | x_{0} - γ (0) |_{1},

i.e. the solution converges to γ with exponential rate α.

Complete proofs of theorems A.1 and A.2 are provided in the following seven subsections. We begin by showing in lemma A.4 that solutions of $\dot{x} = f (t, x)$ , $f \in F$ that start in the closed state-space $R_{+}^{N}$ are unique and remain in $R_{+}^{N}$ for all positive times. In the second subsection, we prove that for the subset of linear vector fields f in $F$ the flow is cooperative, and non-expansive or even contractive in the case of irreducibility. The latter property is then generalized to the nonlinear setting (theorem A.8), which is enough to prove theorem A.2 in subsection A.4. The cooperative behaviour for nonlinear vector fields is stated in theorem A.9. In subsection A.6, we argue that the non-expansiveness of the flow together with the existence of a fixed point in the ω-limit set of the period map implies the asymptotic periodicity of the solution. The proof that such a fixed point exists is deferred to the final subsection. It uses the cooperative behaviour of the flow as well as the fact that the first integral H has a positive gradient, i.e. ∇H∈ int $(R_{+}^{N})$ .

A.1. Positive invariance of $R_{+}^{N}$

Our first goal is to establish in lemma A.4 below that for any $x_{0} \in R_{+}^{N}$ a unique solution x(t;x₀) exists for all $t \in [0, \infty)$ and remains in the closed cone $R_{+}^{N}$ . Denote

B := {B \in R^{N \times N} | B is Metzler and \sum_{i = 1}^{N} B_{i j} = 0 for all j \in {1, \dots, N}} .

Proposition A.3

Assume that $f \in F$ . Then there exists a continuous map $B : [0, \infty) \times R_{+}^{N} \times R_{+}^{N} \to B$ such that

f (t, x) - f (t, y) = B (t, x, y) (x - y) for all t \geq 0 and all x, y \in R_{+}^{N} .

A 2

Moreover, for any indexjthe following property holds. If

x \in R_{+}^{N}

withx_j=0 thenf_j(t,x)≥0 for allt≥0.

Proof.

Equation (A 2) follows from the fundamental theorem of calculus with

B (t, x, y) = \int_{0}^{1} J (t, y + s (x - y)) d s for all t \geq 0 and x, y \in int (R_{+}^{N}) .

By assumptions (iii) and (iv) on f, we conclude that

B (t, x, y) \in B

. Moreover, assumption (ii) implies that this formula actually defines B as a continuous map on

[0, \infty) \times R_{+}^{N} \times R_{+}^{N}

into

B

. Equation (A 2) also holds on the closed cone

R_{+}^{N}

due to assumption (i). The second claim follows from f(t,x)=B(t,x,0)x and the fact that B(t,x,0) is Metzler. ▪

Lemma A.4

Assume that $f \in F$ . Then for every $x_{0} \in R_{+}^{N}$ the initial value problem $\dot{x} = f (t, x)$ , x(0)=x₀, admits a unique solution $x (\cdot; x_{0}) : [0, \infty) \to R_{+}^{N}$ . Moreover,

H (x) = \sum_{j = 1}^{N} x_{j}

is a first integral of the dynamics, i.e. (d/dt)H(x(t;x₀))≡0.

Proof.

By assumption (i) on f and by proposition A.3, the vector field f satisfies the hypotheses of the Picard–Lindelöf theorem, however, with a domain of definition $R_{+}^{n}$ which is not open. Introduce the auxiliary extension

\tilde{f} (t, x) := f (t, | x_{1} |, \dots, | x_{N} |) .

Note that

\tilde{f} : [0, \infty) \times R^{N} \to R^{N}

is well defined. The Picard–Lindelöf theorem yields existence and uniqueness of a maximal solution

\tilde{x} (\cdot; x_{0}) : J_{x_{0}} \to R^{N}

of the initial value problem

\dot{x} = \tilde{f} (t, x)

, x(0)=x₀, with J_x₀=[0,b_x₀) for some b_x₀>0 that might be infinite. Note that H is a first integral by property (iv) of f, that carries over to

\tilde{f}

On $R_{+}^{N}$ , the solutions x( ⋅ ;x₀) and $\tilde{x} (\cdot; x_{0})$ coincide. We now show that for $x_{0} \in R_{+}^{N}$ the solution $\tilde{x} (t; x_{0}) \in R_{+}^{N}$ for all t∈J_x₀. If x₀=0 then $\tilde{x} (t; x_{0}) \equiv 0$ by assumption (v) on f, so $\tilde{x} (t; x_{0}) \in R_{+}^{N}$ . If $x_{0} \in R_{+}^{N} ∖ {0}$ , we argue by contradiction. Assume that there exists τ>0 such that $\tilde{x} (τ; x_{0}) \in (R^{N} ∖ R_{+}^{N})$ . For $ε \in R$ , let

{\tilde{f}}_{ε} : [0, \infty) \times R^{N} \to R^{N}, {\tilde{f}}_{ε} (t, x) := \tilde{f} (t, x) + ε (\begin{matrix} 1 & \dots & 1 \\ ⋮ & ⋮ & ⋮ \\ 1 & \dots & 1 \end{matrix}) x - ε N x

and denote by

\tilde{x} (\cdot; x_{0}, ε)

the unique maximal solution of

\dot{x} = {\tilde{f}}_{ε} (t, x)

, x(0)=x₀. Note that by the definition of

{\tilde{f}}_{ε}

the function H is also a first integral for this dynamical system.

Now by the continuous dependence of solutions on parameters there exists ε₀>0 such that $\tilde{x} (τ; x_{0}, ε_{0}) \in (R^{N} ∖ R_{+}^{N})$ . This implies that there exists a time s∈[0,τ) such that x(s;x₀,ε₀) leaves $R_{+}^{N}$ , i.e. there exists k∈{1,…,N} with

{\tilde{x}}_{k} (s; x_{0}, ε_{0}) = 0 and {\dot{\tilde{x}}}_{k} (s; x_{0}, ε_{0}) \leq 0.

A 3

We will show that this is a contradiction. As

{\dot{\tilde{x}}}_{k} (s; x_{0}, ε_{0}) = {({\tilde{f}}_{ε} (s, \tilde{x} (s; x_{0}, ε_{0})))}_{k}

, the definition of

{\tilde{f}}_{ε}

and (A 3) yield

{\dot{\tilde{x}}}_{k} (s; x_{0}, ε_{0}) = {(\tilde{f} (s, \tilde{x} (s; x_{0}, ε_{0})))}_{k} + ε_{0} \sum_{i = 1}^{N} {\tilde{x}}_{i} (s; x_{0}, ε_{0})

. By the second statement in proposition A.3,

{(\tilde{f} (s, \tilde{x} (s; x_{0}, ε_{0})))}_{k} \geq 0

, so

{\dot{\tilde{x}}}_{k} (s; x_{0}, ε_{0}) \geq ε_{0} \sum_{i = 1}^{N} {\tilde{x}}_{i} (s; x_{0}, ε_{0})

. As H is also a first integral for

\dot{y} = {\tilde{f}}_{ε} (t, y)

, we obtain

{\dot{\tilde{x}}}_{k} (s; x_{0}, ε_{0}) \geq ε_{0} H (x_{0}) > 0

. This contradicts (A 3).

We have established that $\tilde{x}$ remains in the compact set ${z \in R_{+}^{N} | H (z) = H (x_{0})}$ and it follows that the solution $\tilde{x} (t; x_{0})$ exists for all t≥0. As the solutions x( ⋅ ;x₀) and $\tilde{x} (\cdot; x_{0})$ coincide on $R_{+}^{N}$ , this completes the proof. ▪

A.2. Linear time-varying systems

The properties that are essential in the proofs of our main results are cooperativeness, non-expansiveness and contractivity of the flow. As it turns out, it is convenient to first prove these properties for linear time-varying systems. Let

A := {A : [0, \infty) \to B | A is continuous} .

Lemma A.5

Assume that $A \in A$ . Then the initial value problem $\dot{x} = A (t) x$ , $x (0) = x_{0} \in R_{+}^{N}$ , has a unique solution $x (\cdot; x_{0}) : [0, \infty) \to R^{N}$ that satisfies the following properties:

(a) $x (t; x_{0}) \in R_{+}^{N}$ andH(x(t;x₀))=H(x₀) for allt≥0.
(b) Ifx_j(t*;x₀)>0 for somej∈{1,…,N} andt*≥0, thenx_j(t;x₀)>0 for allt≥t*.
(c) Ifx₀≠0 andA(t*) is irreducible for somet*≥0, then $x (t; x_{0}) \in int (R_{+}^{N})$ for allt>t*.

Proof.

Existence and uniqueness of the solution are immediate from the linearity and continuity of A.

The proof of (a) follows from lemma A.4 and the fact that f(t,x):=A(t)x belongs to $F$ for each $A \in A$ .

To prove (b), assume that x_j(t*;x₀)>0. Let y(t):=x_j(t;x₀). Then y solves the scalar initial value problem

\dot{y} = g (t, y), y (t^{*}) = x_{j} (t^{*}; x_{0}) > 0,

where

g (t, y) := a_{j j} (t) y + b (t), b (t) := \sum_{k \neq j} a_{j k} (t) x_{k} (t; x_{0}) \geq 0.

Thus, letting

q (t) := \int_{t^{*}}^{t} a_{j j} (u) d u

yields

\begin{aligned} y (t) & = e^{q (t)} (y (t^{*}) + \int_{t^{*}}^{t} e^{- q (s)} b (s) d s) \\ \geq e^{q (t)} y (t^{*}) \\ > 0, \end{aligned}

for all t≥t*.

To prove property (c), first note that irreducibility of A(t*) implies that there exists δ>0 such that

A (t) is irreducible for all t \in [t^{*}, t^{*} + δ) .

A 4

This follows from the fact that irreducibility is equivalent to the associated adjacency graph being strongly connected, i.e. certain edges have positive weights, and the continuity of A(t).

Pick $x_{0} \in R_{+}^{n} ∖ {0}$ . We consider two cases.

Case 1: $x (t^{*}; x_{0}) \in int (R_{+}^{N})$ . Then the claim follows from property (b).

Case 2: $x (t^{*}; x_{0}) \in \partial R_{+}^{N}$ . Fix t>t*. As H(x(t*;x₀))=H(x₀)>0, there exists k∈{1,…,N−1} such that exactly k entries of x(t*;x₀) are positive and the other entries are zero. Note that by property (b), these k entries remain positive for all t≥t*. Assume w.l.o.g. that the first k entries of x(t*;x₀) are positive. Then in block form

\dot{x} (t^{*}; x_{0}) = (\begin{array}{cc} U & V \\ Y & Z \end{array}) (\begin{matrix} x_{1} (t^{*}; x_{0}) \\ ⋮ \\ x_{k} (t^{*}; x_{0}) \\ 0 \end{matrix}),

with

U \in R^{k \times k}

and

Z \in R^{(N - k) \times (N - k)}

. As A(t*) is Metzler and irreducible, every entry of Y is non-negative and at least one entry is positive, so there exists j>k such that

{\dot{x}}_{j} (t^{*}; x_{0}) > 0

. Therefore, at least k+1 entries of x(t;x₀) are positive for t>t*. Now an inductive argument and using (A 4) completes the proof. ▪

Let

S := {Q \in R^{N \times N} | q_{j k} \geq 0 for all j, k = 1, \dots, N and \sum_{j = 1}^{N} q_{j k} = 1 for all k = 1, \dots, N}

denote the set of N×N stochastic matrices, and let

S^{+} := {Q \in S | q_{j k} > 0 for all j, k = 1, \dots, N}

denote the subset of stochastic matrices with positive entries. For

A \in A

let

Φ_{A} : [0, \infty) \to R^{N \times N}

be the fundamental matrix of

\dot{x} = A x

, that is, the solution of

\dot{Φ} (t) = A (t) Φ (t), Φ (0) = I_{N} .

As the columns of Φ_A(t) are x(t;e^j), where

e^{j} \in R_{+}^{N}

denotes the jth canonical unit vector, the next result follows from properties (a) and (c) in lemma A.5.

Corollary A.6

Assume that $A \in A$ . Then

(a) $Φ_{A} (t) \in S$ for allt≥0.
(b) IfA(t*) is irreducible for somet*≥0, then $Φ_{A} (t) \in S^{+}$ for allt>t*.

We now use this to prove non-expansiveness with respect to the ℓ₁-norm and contractivity in the case of irreducibility. The first step is to note that stochastic matrices have useful properties with respect to this norm.

Proposition A.7

If $Q \in S$ and $x \in R^{N}$ , then |Qx|₁≤|x|₁. If $Q \in S^{+}$ and $x \in (R^{n} ∖ {0})$ withH(x)=0, then |Qx|₁<|x|₁.

Proof.

The first statement follows from

\begin{aligned} | Q x |_{1} & = \sum_{j = 1}^{N} | \sum_{k = 1}^{N} Q_{j k} x_{k} | \\ \leq \sum_{j = 1}^{N} \sum_{k = 1}^{N} Q_{j k} | x_{k} | \\ = \sum_{k = 1}^{N} | x_{k} | \\ = | x |_{1} . \end{aligned}

A 5

To prove the second statement, pick x≠0 such that H(x)=0. Then there exist k₁,k₂∈{1,…,N} such that x_k₁<0<x_k₂. Thus, if $Q \in S^{+}$ then for any j∈{1,…,N},

| \sum_{k = 1}^{N} Q_{j k} x_{k} | < \sum_{k = 1}^{N} Q_{j k} | x_{k} |

because the sum on the left contains both positive and negative terms. Now arguing as in (A 5) completes the proof. ▪

A.3. Non-expansiveness and contractivity

Using the results for time-varying linear systems, we now turn to proving non-expansiveness and contractivity for the nonlinear dynamical system.

Theorem A.8

Suppose that $f \in F$ . Then

(a) For any $x_{0}, y_{0} \in R_{+}^{N}$ the function
$t \mapsto | x (t; x_{0}) - x (t; y_{0}) |_{1} is non-increasing on [0, \infty) .$ A 6
(b) If there exists t*≥0 such that J(t*,x) is irreducible for all x∈Ω then for any $\hat{t} > t^{*}$ there exists $M_{\hat{t}} < 1$ such that
$| x (\hat{t}; x_{0}) - x (\hat{t}; y_{0}) |_{1} \leq M_{\hat{t}} | x_{0} - y_{0} |_{1} for all x_{0}, y_{0} \in Ω .$ A 7

Proof.

Fix t₁≥0 and x₀,y₀∈Ω. Let z(t):=x(t+t₁;x₀)−x(t+t₁;y₀). By proposition A.3, $\dot{z} (t) = C (t) z (t)$ , with z(0)= z₀:=x(t₁;x₀)−x(t₁;y₀), and $C \in A$ given by

C (t) := B (t + t_{1}, x (t + t_{1}; x_{0}), x (t + t_{1}; y_{0})),

for all t≥0. Hence, z(t)=Φ_C(t)z₀ for all t≥0.

To prove (A 6), pick t₂≥t₁. By corollary A.6(a), $Φ_{C} (t_{2} - t_{1}) \in S$ and proposition A.7 yields

\begin{aligned} | x (t_{2}; x_{0}) - x (t_{2}; y_{0}) |_{1} & = | z (t_{2} - t_{1}) |_{1} \\ = | Φ_{C} (t_{2} - t_{1}) z_{0} |_{1} \\ \leq | z_{0} |_{1} \\ = | x (t_{1}; x_{0}) - x (t_{1}; y_{0}) |_{1} \end{aligned}

and this proves (A 6).

To prove (A 7), let

P := {x \in R^{N} | H (x) = 0 and | x |_{1} = 1} .

Note that if y,z∈Ω with y≠z then

(y - z) / | y - z |_{1} \in P

Pick $\hat{t} > t^{*}$ and x₀,y₀∈Ω. Equation (A 7) clearly holds if x₀=y₀, so we may assume that x₀≠y₀. Then

\begin{aligned} | x (\hat{t}; x_{0}) - x (\hat{t}; y_{0}) |_{1} & = | Φ_{C} (\hat{t}) (x_{0} - y_{0}) |_{1} \\ = m (x_{0}, y_{0}, \frac{x_{0} - y_{0}}{| x_{0} - y_{0} |_{1}}) | x_{0} - y_{0} |_{1}, \end{aligned}

where

m : Ω \times Ω \times P \to R

is defined by

m (x_{0}, y_{0}, v) := | Φ_{C} (\hat{t}) v |_{1}

. Let

M_{\hat{t}} := max {m (x, y, v) | x \in Ω, y \in Ω, v \in P} .

This is well defined, as the matrix C(t)=C_x₀,y₀(t)=B(t,x(t;x₀),x(t;y₀)) depends continuously on x₀,y₀∈Ω for all t≥0 (see proposition A.3) and thus

Φ_{C_{x_{0}, y_{0}}} (\hat{t})

is also continuous in x₀,y₀. We conclude that

| x (\hat{t}; x_{0}) - x (\hat{t}; y_{0}) |_{1} \leq M_{\hat{t}} | x_{0} - y_{0} |_{1} .

Thus, to complete the proof we only need to show that

M_{\hat{t}} < 1

. To prove this, denote x*:=x(t*;x₀), y*:=x(t*;y₀). Then

C_{x_{0}, y_{0}} (t^{*}) = B (t^{*}, x^{*}, y^{*}) = \int_{0}^{1} J (t^{*}, y^{*} + s (x^{*} - y^{*})) d s .

As J is Metzler (i.e. all its off-diagonal elements are non-negative) and, by assumption, J(t*,z) is irreducible for all z∈Ω, we conclude that C_x₀,y₀(t*) is also irreducible. Corollary A.6 implies that

Φ_{C_{x_{0}, y_{0}}} (\hat{t}) \in S^{+}

. Picking a maximizer

(x_{0}, y_{0}, v_{0}) \in Ω \times Ω \times P

of m, i.e.

M_{\hat{t}} = m (x_{0}, y_{0}, v_{0})

, it follows from proposition A.7 that

M_{\hat{t}} < | v_{0} |_{1} = 1

. ▪

A.4. Proof of theorem A.2

We can now prove theorem A.2. We note that the proof proceeds without the explicit use of the cooperative behaviour of dynamical systems (though we will use this property for the proof of theorem A.1; see subsections A.5 and A.6).

Note that for $f \in F_{T}$ a solution $\dot{γ} = f (t, γ)$ , $γ : [0, \infty) \to Ω$ , is T-periodic if and only if x(T;γ(0))=γ(0). Thus, consider the period map $P_{T} : R_{+}^{N} \to R_{+}^{N}$ defined by

P_{T} (a) := x (T; a) .

A 8

In other words, P_T(a) is the value of x(T) for the initial condition x(0)=a. Observe that P_T(Ω)⊆Ω (as H is a first integral of

\dot{x} = f (t, x)

). Moreover, for

f \in (F_{T} \cap F_{irr}^{Ω})

there exists t*∈[0,T) such that J(t*,x) is irreducible for all x∈Ω. Then T>t*, so theorem A.8(b) implies that P_T is Lipschitz on the closed set Ω with Lipschitz constant M_T<1. The Banach fixed point theorem implies that P_T has a unique fixed point in Ω, that is, there exists a unique T-periodic function

γ : R \to Ω

that solves

\dot{γ} = f (t, γ)

. Fix α>0 such that

max {\frac{1}{2}, M_{T}} \leq e^{- α T} .

A 9

Pick x₀∈Ω and t≥0, and let

k \in N_{0}

be such that kT≤t≤(k+1)T. Then theorem A.8 yields

\begin{aligned} | x (t; x_{0}) - γ (t) |_{1} & \leq | x (k T; x_{0}) - γ (k T) |_{1} \\ \leq {(M_{T})}^{k} | x_{0} - γ (0) |_{1} \\ \leq e^{- α T k} | x_{0} - γ (0) |_{1} \\ = e^{- α T (k + 1)} e^{α T} | x_{0} - γ (0) |_{1} \\ \leq e^{α T} e^{- α t} | x_{0} - γ (0) |_{1} \\ \leq 2 e^{- α t} | x_{0} - γ (0) |_{1}, \end{aligned}

where the last inequality follows from (A 9). Thus, every solution x(t;x₀) converges to the unique periodic solution γ(t) at an exponential rate.

A.5. Cooperative behaviour

For the proof of theorem A.1, we use some elegant topological ideas from [44,45], which are based on the cooperative behavior of the dynamical system generated by a vector field $f \in F$ . To formulate this concept, we introduce some more notation.

Let $p, q \in R^{N}$ , and let A be a bounded non-empty subset of $R^{N}$ . Then we write

(a) p≤q :⇔ p_j≤q_j for all j∈{1,…,N},
(b) $[p, q] := {x \in R^{N} | p \leq x \leq q}$ ,
(c) $inf A := c \in R^{N}$ with $c_{j} := inf {a_{j} | a \in A}$ , $sup A := d \in R^{N}$
with $d_{j} := sup {a_{j} | a \in A}$ ,
(d) p≤A ⇔ p≤a for all a∈A
q≥A ⇔ q≥a for all a∈A.

It is straightforward to verify that $inf A \leq A \leq sup A$ and that for all $x, y \in R^{N}$ with x≤A, y≥A we have $x \leq inf A$ and $y \geq sup A$ .

The following theorem summarizes the monotone behaviour with respect to the order ≤.

Theorem A.9

Let $f \in F$ and $x_{0}, y_{0} \in R_{+}^{N}$ with x₀≤y₀. Then x(t;x₀)≤x(t;y₀) for all t≥0. If, in addition, (x₀)_j<(y₀)_jfor some j∈{1,…,N}, then x_j(t;x₀)<x_j(t;y₀) for all t≥0.

Proof.

We use again that z(t):=x(t;y₀)−x(t;x₀) satisfies $\dot{z} = C (t) z$ , $z (0) = z_{0} := y_{0} - x_{0} \in R_{+}^{N}$ with $C (t) := B (t, x (t; y_{0}), x (t; x_{0})) \in B$ and $C \in A$ . The claim then follows from statements (a) and (b) of lemma A.5. ▪

A.6. ω-Limit sets of P_T and proof of theorem A.1

The concept of ω-limit sets for the discrete time dynamical system induced by the period map P_T:Ω→Ω from (A 8) is pivotal for the proof of theorem A.1. For a∈Ω this set is

ω_{T} (a) := {x \in Ω | there is a sequence n_{k} \to \infty with lim_{k \to \infty} P_{T}^{n_{k}} (a) = x} .

We first state a few standard facts about this set.

Proposition A.10

Let $f \in F$ . PickT>0. Then for alla∈Ωwe have

(a) ω_T(a) is a closed, non-empty subset ofΩ;
(b) P_T(ω_T(a))=ω_T(a);
(c) ω_T(b)⊆ω_T(a) for allb∈ω_T(a).

Proof.

Statements (a) and (b) follow from [46], eqn. (4.7.2) and lemma 4.7.4 because $P_{T}^{k} (a)$ evolves in the compact set Ω.

To prove (c), pick b∈ω_T(a). Property (b) implies that $P_{T}^{n_{k}} (b) \in ω_{T} (a)$ for all $n_{k} \in N$ , and as ω_T(a) is closed $lim_{k \to \infty} P_{T}^{n_{k}} (b) \in ω_{T} (a)$ . ▪

The following lemma, that is proved in the subsequent subsection, provides all that is needed to prove theorem A.1.

Lemma A.11

Let $f \in F_{T}$ . Then for everyx∈Ωits limit setω_T(x) contains a fixed point ofP_T.

Proof of theorem A.2

Pick x₀∈Ω and denote by z∈ω_T(x₀) the fixed point of P_T that exists according to lemma A.11. As z is a fixed point and as $f \in F_{T}$ is T-periodic, the solution γ(t):=x(t;z) is also T-periodic. As z∈ω_T(x₀) there exists a subsequence $P_{T}^{n_{k}} (x_{0}) \to z$ as $k \to \infty$ . For any $j \in N$ and t≥n_jT, we have by theorem A.8(a),

| x (t; x_{0}) - γ (t) |_{1} \leq | x (T n_{j}; x_{0}) - γ (T n_{j}) |_{1} = | P_{T}^{n_{j}} (x_{0}) - z |_{1} .

t \to \infty

, we can take

j \to \infty

and this yields

lim_{t \to \infty} | x (t; x_{0}) - γ (t) |_{1} = 0,

proving theorem A.1. ▪

Note that the proof of theorem A.1 shows that for any x₀∈Ω the ω-limit set ω_T(x₀) cannot contain more than one fixed point of P_T. Thus, the statement in lemma A.11 can actually be strengthened to ω_T(x₀) contains exactly one fixed point of P_T.

The next subsection contains the proof of the crucial lemma A.11. We have adapted the ideas presented by Ji-Fa in [45] to our setting, which led to somewhat simplified arguments. In particular, we can replace [44], proposition 1 that is used in the proof of lemma 3.2 of [45] by a standard application of Brouwer’s fixed point theorem. Indeed, the claim of lemma A.11 is the existence of a fixed point of the map P_T in ω_T(x). We show that ω_T(x) contains an element y so that ω_T(y) is a singleton, say, ω_T(y)={z}. Then z∈ω_T(x) by proposition A.10(c) and z is a fixed point by proposition A.10(b).

A.7. Proof of lemma A.11

We use the following notation. For y∈Ω, let

\begin{aligned} p (y) & := inf ω_{T} (y), \\ q (y) & := sup ω_{T} (y) . \end{aligned}

Note that as

ω_{T} (y) \subseteq Ω \subset R_{+}^{N}

, p(y) and q(y) exist in

R_{+}^{N}

. Moreover, p(y)≤q(y) and ω_T(y) is a singleton if and only if p(y)=q(y). The existence of the desired element y in ω_T(x) is established by contradiction. Suppose that no such y exists. Denote by y₀ an element in ω_T(x) that minimizes the number of coordinates j for which (p(y))_j and (q(y))_j differ. We then show that there exists z₀∈ω_T(x) for which p(z₀) and q(z₀) differ in a smaller number of coordinates than p(y₀) and q(y₀). Part (c) of lemma A.12 below states a fact that is essential for the construction of z₀. What is also crucial for the proof is the observation that p(y) and q(y) are fixed points of P_T which is formulated in lemma A.12(a). Its short proof demonstrates why monotone dynamical systems admitting a first integral with a positive gradient are special.

For y∈Ω, let

J_{y} := {j \in {1, \dots, N} | {(p (y))}_{j} \neq {(q (y))}_{j}},

A 10

i.e. the set of indices for which (p(y))_j and (q(y))_j differ.

Lemma A.12

Let $f \in F_{T}$ . Then for anyy∈Ω,

(a) p(y), q(y) are fixed points ofP_T.
(b) P_T([p(y),q(y)])⊆[p(y),q(y)].
(c) Ifp(y)≠q(y) then for anyz∈ω_T(y) there exists aj∈J_ysuch thatz_j=(p(y))_j.

To explain property (c), we introduce more notation. For $v, w \in R^{N}$ let

Δ (v, w) := # {j \in {1, \dots, N} | v_{j} \neq w_{j}} .

Consider the case p(y)≠q(y). For any j∈{1,…,N}∖J_y, we have p_j(y)=q_j(y), so z_j=p_j=q_j. For any j∈J_y, we have p_j(y)<q_j(y) and property (c) implies that there exists at least one such index such that p_j(y)=z_j. We conclude that

Δ (z, p (y)) < Δ (q (y), p (y)) .

A 11

Proof of lemma A.12

To simplify the notation, we write from hereon p,q for p(y),q(y). Pick z∈ω_T(y). As p≤z, theorem A.9 implies that P_T(p)≤P_T(z). Thus, P_T(p)≤ω_T(y) by proposition A.10(b), and consequently P_T(p)≤p. As H(P_T(p))=H(p), it follows that P_T(p)=p. The proof that P_T(q)=q proceeds analogously. This proves claim (a).

Claim (b) follows directly from (a) and theorem A.9.

We prove (c) by contradiction. Assume that there exists z∈ω_T(y) such that z_j>p_j for all j∈J_y. Set

ε := min {z_{j} - p_{j} | j \in J_{y}} .

A 12

Then ε>0. Define

Γ := {x \in [p, q] | H (x) = H (p) + \frac{ε}{2}} .

A 13

Note that for any j∈{1,…,N}∖J_y we have p_j=q_j so if x∈Γ then x_j=p_j=q_j. The set Γ is not empty, as H(q)≥H(z)≥H(p)+ε. Γ is also compact and convex. Statement (b) and the fact that H is a first integral imply that P_T(Γ)⊆Γ. The Brouwer fixed point theorem thus yields the existence of a fixed point of P_T in Γ, that is, there exists x*∈Γ such that P_T(x*)=x*.

As z∈ω_T(y) there exists $n^{*} \in N$ such that $| P_{T}^{n^{*}} (y) - z |_{1} < ε / 3 N$ . Define y* by

y_{j}^{*} := {\begin{cases} {(P_{T}^{n^{*}} (y))}_{j}, & j \in J_{y} \\ z_{j} = p_{j} = q_{j}, & otherwise . \end{cases}

| P_{T}^{n^{*}} (y) - y^{*} |_{1} + | y^{*} - z |_{1} \leq | P_{T}^{n^{*}} (y) - z |_{1}

, we have

| P_{T}^{n^{*}} (y) - y^{*} |_{1} < ε / 3 N

and

| y^{*} - z |_{1} < \frac{ε}{3 N} .

A 14

Observe that

x_{j}^{*} = p_{j} = z_{j} = y_{j}^{*} for j \in {1, \dots, N} ∖ J_{y}

and

x_{j}^{*} \leq p_{j} + \frac{ε}{2} \leq z_{j} - \frac{ε}{2} < y_{j}^{*} for j \in J_{y},

where the first inequality follows from (A 13), the second from (A 12) and the third from (A 14).

Summarizing, x*≤y*. As P_T(x*)=x*, theorem A.9 implies that $x^{*} \leq P_{T}^{k} (y^{*})$ for all k≥0. As $\sum_{j} (x_{j}^{*} - p_{j}) = ε / 2$ , there exists j₀∈{1,…,N} such that $x_{j_{0}}^{*} \geq p_{j_{0}} + ε / 2 N$ . Then ${(P_{T}^{k} (y^{*}))}_{j_{0}} \geq p_{j_{0}} + ε / 2 N$ for all k≥0.

As P_T is non-expansive by theorem A.8(a), we have in addition $| P_{T}^{k} (y^{*}) - P_{T}^{k + n^{*}} (y) |_{1} \leq | y^{*} - P_{T}^{n^{*}} (y) |_{1} < ε / 3 N$ for all k≥0. Hence for all m≥n*,

{(P_{T}^{m} (y))}_{j_{0}} \geq p_{j_{0}} + \frac{ε}{2 N} - \frac{ε}{3 N} = p_{j_{0}} + \frac{ε}{6 N} .

We conclude that any r∈ω_T(y) satisfies r_j₀≥p_j₀+ε/6N, and thus

p_{j_{0}} = inf {r_{j_{0}} | r \in ω_{T} (y)} \geq p_{j_{0}} + \frac{ε}{6 N}

and this contradiction proves (c). ▪

We can now prove the crucial lemma A.11.

Proof of lemma A.11

As noted above, we need to show that there exists y∈ω_T(x) such that p(y)=q(y). To do this, for y∈Ω let

m (y) := Δ (p (y), q (y))

and for x∈Ω, let

α (x) := min {m (y) | y \in ω_{T} (x)} .

It suffices to show that α(x)=0 for all x∈Ω. We achieve this by contradiction. Assume that there exists x∈Ω for which α(x)>0. Then there exists y₀∈ω_T(x) with α(x)=m(y₀)>0. Let

\tilde{M} := max {Δ (z, p (y_{0})) | z \in ω_{T} (y_{0})}

. Then (A 11) yields

\tilde{M} < α (x) .

A 15

Choose z₀∈ω_T(y₀) such that

Δ (z_{0}, p (y_{0})) = \tilde{M}

. We now show that

m (z_{0}) \leq \tilde{M}

. To this end, define

J := {j \in {1, \dots . N} | {(z_{0})}_{j} > p_{j} (y_{0})} .

Note that

\tilde{M} = # J

. As z₀∈[p(y₀),q(y₀)], we have J⊆J_y₀ and by lemma A.12 we have J≠J_y₀. Consider the sequence

z (k) := P_{T}^{k} (z_{0}) \in ω_{T} (y_{0})

. By the second statement of theorem A.9, we have z_j(k)>p_j(y₀) for all j∈J and all

k \in N

(recall that p(y₀) is a fixed point). Hence

Δ (z (k), p (y_{0})) \geq # J = \tilde{M}

. As z(k)∈ω_T(y₀), we have by the definition of

\tilde{M}

also

Δ (z (k), p (y_{0})) \leq \tilde{M}

and, therefore,

Δ (z (k), p (y_{0})) = \tilde{M}

, and

J = {j \in {1, \dots, N} | z_{j} (k) > p_{j} (y_{0})} for all k \in N .

Thus, z_j(k)=p_j(y₀) for all j∈{1,…,N}∖J, so

ω_{T} (z_{0}) \subseteq {v \in R^{N} | v_{j} = p_{j} (y_{0}) for all j \in {1, \dots, N} ∖ J} .

This implies that p_j(z₀)=q_j(z₀) for all j∈{1,…,N}∖J, and thus

m (z_{0}) = Δ (p (z_{0}), q (z_{0})) \leq # J = \tilde{M}

. Combining this with the fact that z₀∈ω_T(y₀)⊆ω_T(x) and (A 15) yields

α (x) \leq m (z_{0}) \leq \tilde{M} < α (x)

and this contradiction completes the proof of lemma A.11. ▪

Footnotes

1 That is, the real part of every eigenvalue of A is negative.

2 We consider the continuous time version of TASEP here.

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.

References

1
Gillespie DT. 1992A rigorous derivation of the chemical master equation. Phys. A: Stat. Mech. Appl. 188, 404–425. (doi:10.1016/0378-4371(92)90283-V) Crossref, ISI, Google Scholar
2
Haag G. 2017Modelling with the master equation: solution methods and applications in social and natural sciences. Cham, Switzerland: Springer International Publishing. Crossref, Google Scholar
3
Van Kampen NG. 2007Stochastic processes in physics and chemistry. 3rd edn. Amsterdam, The Netherlands: Elsevier. Google Scholar
4
Resnick SI. 2002Adventures in stochastic processes. Boston, MA: Birkhauser. Crossref, Google Scholar
5
Toral R, Colet P. 2014Stochastic numerical methods. Weinheim, Germany: Wiley. Crossref, Google Scholar
6
Schadschneider A, Chowdhury D, Nishinari K. 2011Stochastic transport in complex systems: from molecules to vehicles. Amsterdam, The Netherlands: Elsevier. Google Scholar
7
Krug J. 2016Nonequilibrium stationary states as products of matrices. J. Phys. A: Math. Theor. 49, 421002. (doi:10.1088/1751-8113/49/42/421002) Crossref, ISI, Google Scholar
8
Nadler W, Schulten K. 1986Generalized moment expansion for observables of stochastic processes in dimensions d>1: application to Mossbauer spectra of proteins. J. Chem. Phys. 84, 4015–4025. (doi:10.1063/1.450061) Crossref, ISI, Google Scholar
9
Russo G, di Bernardo M, Sontag ED. 2010Global entrainment of transcriptional systems to periodic inputs. PLoS Comput. Biol. 6, e1000739. (doi:10.1371/journal.pcbi.1000739) Crossref, PubMed, ISI, Google Scholar
10
Keith WL, Rand RH. 19841:1 and 2:1 phase entrainment in a system of two coupled limit cycle oscillators. J. Math. Biol. 20, 133–152. (doi:10.1007/BF00285342) Crossref, ISI, Google Scholar
11
Grassly NC, Fraser C. 2006Seasonal infectious disease epidemiology. Proc. R. Soc. B: Biol. Sci. 273, 2541–2550. (doi:10.1098/rspb.2006.3604) Link, ISI, Google Scholar
12
Donner R. 2009Emergence of synchronization in transportation networks with biologically inspired decentralized control. In Recent advances in nonlinear dynamics and synchronization (eds K Kyamakya, H Unger, JC Chedjou, NF Rulkov, Z Li), pp. 237–275. Studies in Computational Intelligence, vol. 254. Berlin, Germany: Springer-Verlag. Google Scholar
13
Groß D, Arghir C, Dörfler F. 2016On the steady-state behavior of a nonlinear power system model. (http://arxiv.org/abs/1607.01575) Google Scholar
14
Nikolaev EV, Rahi SJ, Sontag ED. 2018Subharmonics and chaos in simple periodically forced biomolecular models. Biophys. J. 114, 1232–1240. (doi:10.1016/j.bpj.2018.01.006) Crossref, PubMed, ISI, Google Scholar
15
MacDonald CT, Gibbs JH, Pipkin AC. 1968Kinetics of biopolymerization on nucleic acid templates. Biopolymers 6, 1–25. (doi:10.1002/bip.1968.360060102) Crossref, PubMed, ISI, Google Scholar
16
Zia R, Dong J, Schmittmann B. 2011Modeling translation in protein synthesis with TASEP: a tutorial and recent developments. J. Stat. Phys. 144, 405–428. (doi:10.1007/s10955-011-0183-1) Crossref, ISI, Google Scholar
17
Shaw LB, Zia RKP, Lee KH. 2003Totally asymmetric exclusion process with extended objects: a model for protein synthesis. Phys. Rev. E 68, 021910. (doi:10.1103/PhysRevE.68.021910) Crossref, PubMed, ISI, Google Scholar
18
Kriecherbauer T, Krug J. 2010A pedestrian’s view on interacting particle systems, KPZ universality, and random matrices. J. Phys. A: Math. Theor. 43, 403001. (doi:10.1088/1751-8113/43/40/403001) Crossref, ISI, Google Scholar
19
Chowdhury D, Santen L, Schadschneider A. 1999Vehicular traffic: a system of interacting particles driven far from equilibrium. Curr. Sci. 77, 411–419. ISI, Google Scholar
20
Margaliot M, Sontag ED. 2018Revisiting totally positive differential systems: a tutorial and new results. Submitted. [Online]. See https://arxiv.org/abs/1802.09590. Google Scholar
21
Aminzare Z, Sontag ED. 2014Contraction methods for nonlinear systems: a brief introduction and some open problems. In Proc. 53rd IEEE Conf. on Decision and Control, Los Angeles, CA, pp. 3835–3847. Piscataway, NJ: IEEE. (doi:10.1109/CDC.2014.7039986) Google Scholar
22
Lohmiller W, Slotine J-JE. 1998On contraction analysis for non-linear systems. Automatica 34, 683–696. (doi:10.1016/S0005-1098(98)00019-3) Crossref, ISI, Google Scholar
23
Margaliot M, Coogan S. 2017Approximating the frequency response of contractive systems, ArXiv e-prints. [Online]. See http://adsabs.harvard.edu/abs/2017arXiv170206576M. Google Scholar
24
Zadeh LA, Desoer CA. 1963Linear system theory. New York, NY: McGraw-Hill. Google Scholar
25
Smith HL. 1995Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems. Mathematical Surveys and Monographs, vol. 41. Providence, RI: American Mathematical Society. Google Scholar
26
Mierczynski J. 1991A class of strongly cooperative systems without compactness. Colloq. Math. 62, 43–47. (doi:10.4064/cm-62-1-43-47) Crossref, Google Scholar
27
Margaliot M, Tuller T, Sontag ED. 2017Checkable conditions for contraction after small transients in time and amplitude. In Feedback stabilization of controlled dynamical systems: in honor of Laurent Praly (ed. N Petit), pp. 279–305. Cham, Switzerland: Springer International Publishing. Google Scholar
28
Horn RA, Johnson CR. 2013Matrix analysis, 2nd edn. Cambridge, UK: Cambridge University Press. Google Scholar
29
Popkov V, Salerno M, Schütz GM. 2008Asymmetric simple exclusion process with periodic boundary driving. Phys. Rev. E 78, 011122. (doi:10.1103/PhysRevE.78.011122) Crossref, ISI, Google Scholar
30
Yesil AF, Yalabik MC. 2016Dynamical phase transitions in totally asymmetric simple exclusion processes with two types of particles under periodically driven boundary conditions. Phys. Rev. E 93, 012123. (doi:10.1103/PhysRevE.93.012123) Crossref, PubMed, ISI, Google Scholar
31
Basu U, Chaudhuri D, Mohanty PK. 2011Bimodal response in periodically driven diffusive systems. Phys. Rev. E 83, 031115. (doi:10.1103/PhysRevE.83.031115) Crossref, ISI, Google Scholar
32
Margaliot M, Sontag ED, Tuller T. 2014Entrainment to periodic initiation and transition rates in a computational model for gene translation. PLoS ONE 9, e96039. (doi:10.1371/journal.pone.0096039) Crossref, PubMed, ISI, Google Scholar
33
Poker G, Zarai Y, Margaliot M, Tuller T. 2014Maximizing protein translation rate in the nonhomogeneous ribosome flow model: a convex optimization approach. J. R. Soc. Interface 11, 20140713. (doi:10.1098/rsif.2014.0713) Link, ISI, Google Scholar
34
Raveh A, Margaliot M, Sontag ED, Tuller T. 2016A model for competition for ribosomes in the cell. J. R. Soc. Interface 13, 20151062. (doi:10.1098/rsif.2015.1062) Link, ISI, Google Scholar
35
Poker G, Margaliot M, Tuller T. 2015Sensitivity of mRNA translation. Sci. Rep. 5, 12795. (doi:10.1038/srep12795) Crossref, PubMed, ISI, Google Scholar
36
Zarai Y, Margaliot M, Tuller T. 2016On the ribosomal density that maximizes protein translation rate. PLoS ONE 11, 1–26. (doi:10.1371/journal.pone.0166481) Crossref, ISI, Google Scholar
37
Blythe RA, Evans MR. 2007Nonequilibrium steady states of matrix-product form: a solver’s guide. J. Phys. A: Math. Theor. 40, R333–R441. (doi:10.1088/1751-8113/40/46/R01) Crossref, ISI, Google Scholar
38
Derrida B, Evans MR. 1997The asymmetric exclusion model: exact results through a matrix approach. In Nonequilibrium statistical mechanics in one dimension (ed. V Privman), pp. 277–304. Cambridge, UK: Cambridge University Press. Google Scholar
39
Nåsel I. 2011Extinction and quasi-stationarity in the stochastic logistic SIS model. Lecture Notes in Mathematics, vol. 2022. Berlin, Germany: Springer. Google Scholar
40
Bacaër N. 2015On the stochastic SIS epidemic model in a periodic environment. J. Math. Biol. 71, 491–511. (doi:10.1007/s00285-014-0828-1) Crossref, PubMed, ISI, Google Scholar
41
Soper HE. 1929The interpretation of periodicity in disease prevalence. J. R. Stat. Soc. 92, 34–73. (doi:10.2307/2341437) Crossref, Google Scholar
42
Krug J. 1991Boundary-induced phase transitions in driven diffusive systems. Phys. Rev. Lett. 67, 1882–1885. (doi:10.1103/PhysRevLett.67.1882) Crossref, PubMed, ISI, Google Scholar
43
Nakajima F. 1979Periodic time dependent gross-substitute systems. SIAM J. Appl. Math. 36, 421–427. (doi:10.1137/0136032) Crossref, ISI, Google Scholar
44
Dancer EN, Hess P. 1991Stability of fixed points for order-preserving discrete-time dynamical systems. J. Reine Angew. Math. 419, 125–139. (doi:10.1515/crll.1991.419.125) ISI, Google Scholar
45
Ji-Fa J. 1996Periodic monotone systems with an invariant function. SIAM J. Math. Anal. 27, 1738–1744. (doi:10.1137/S003614109326063X) Crossref, ISI, Google Scholar
46
Michel AN, Wang K, Hu B. 2001Qualitative theory of dynamical systems, 2nd edn. Monographs and Textbooks in Pure and Applied Mathematics, vol. 239. New York, NY: Marcel Dekker. Google Scholar