Agent-based modeling: Population limits and large timescales

Social dynamics and collective social phenomena arise in systems of multiple agents that act and interact, often based on incomplete information, within a social network embedded in a common environment, which can be given by geographical conditions, infrastructure, resources, etc., as well as by norms, rules, and narratives. An agent can represent an individual, a household, firm, political or administrative organization, or any type of discrete entity. A wide range of applications, such as innovation spreading (e.g., Ref. 1) or infection kinetics (e.g., Ref. 2) is addressed, using a wide spectrum of methods from data-based micro-simulation of synthetic populations (e.g., Ref. 3) to abstract individual-⁴ and agent-based models (ABMs) for studying underlying dynamic mechanisms. Comparatively, compact models of social dynamics are present especially in the fields of socio- and econonophysics (see, e.g., the recent special issue “From statistical physics to social sciences”).⁵ A prime application lies in opinion dynamics, a field that goes back to the introduction of the voter model by Clifford and Sudbury⁶ in the 1970s, and later on obtained its name by Holley.⁷ The basic idea is that agents imitate the opinion of neighbors. Various modifications in the representation of opinions, imitation details, and the interaction structures exist (see, e.g., Refs. 8–10 for recent overviews). Throughout the paper, we focus on basic ABMs of this kind.

While the focus on the individuals in social systems and their interactions was proposed as early as 1957 by Orcutt,¹¹ only in the 1990s, motivated by the increasing computing power, computer simulations of the resulting dynamics became a commonly used tool. Also nowadays, many models are implemented with the help of (often object-oriented and discrete-time) software frameworks for agent-based modeling (see, e.g., Ref. 12). Not always the model is precisely specified beforehand. The basics of an ABM are easily communicable, e.g., to stakeholders, as individual (inter)action rules at the microlevel and stepwise simulation of the overall system’s dynamics can be explained without requiring mathematical expertise of the audience. However, the great freedom of the modeler in defining details can make model documentation a challenge.¹³ One point that often remains ambiguous in ABMs with a discrete time step—which is the overwhelming majority—is scheduling, that is, which agent(s) influence(s) which other agent(s) in which order in each time step. As Weimer et al.¹⁴ point out, most textbooks do not deeply reflect upon this topic¹⁵ and models commonly use random interaction orders in each step. In this case, the question arises whether agents “see” the changes already made within the same time step by other agents acting earlier. For parallelization of the ABM simulation, this may lead to problems, however. When considering an ABM to be a representation of a real-world social system, time discretization into steps in which each agent acts once is anyhow somewhat artificial. As an alternative to the discrete-time description of ABMs, there is a variety of approaches using continuous-time Markov processes (see, e.g., Refs. 16 and 17) a corresponding software framework is proposed in Ref. 18. In this paper, we will adopt this continuous-time description, not only because it allows us to circumvent the scheduling problem but mainly because it simplifies the mathematical approach taken herein.

Whether we consider simulation-oriented and rather informal ABMs or more precise mathematical modeling approaches for interacting agent dynamics, we always have to face severe practical and theoretical challenges if the number $N$ of agents and/or the timescale $T$ of interest become large. With increasing $N$ and/or $T$⁠, the computational effort increases drastically (in terms of operations per time step and in terms of the total number of time steps required), which renders simulation-based analysis of the dynamics practically infeasible for ABMs with very large $N$ and/or $T$⁠.

Since the emergence of ABMs, several approaches have been discussed for meeting this challenge. The most prominent ones consider

• the large population limit $N→∞$ for not too large, $N$-independent timescale $T$ by constructing appropriate limit equations that describe the effective dynamics of the ABM but themselves are much less demanding computationally and analytically, or

• the behavior of the ABM on exponentially long timescales $T≍exp⁡(ζN)$ via WKB approximations (after Wentzel, Kramers, and Brillouin) or the large deviation principle (LDP).

There are a lot of results and remarks on the relation between (A) and (B), but the picture is far from being complete. This article contributes to making the picture clearer. In particular, it bridges the gap between (A) and (B) by means of the so-called transfer operator approach. For the sake of comprehensibility, the subsequent considerations will be limited to a voter model like network-based form of ABMs without any spatial component.

The main well-known result with regard to (A), the large population limit, is the insight that, for large $N$⁠, the ABM (described as Markov jump process) can be approximated by a closed system of ordinary differential equations (ODE or mean-field ODE).^19,20 As we will see later, this is a pathwise result, i.e., the rescaled trajectory (e.g., given by opinion percentages) of the ABM converges to the trajectory of the ODE system for $N→∞$ for a fixed finite time interval $[0,T]$⁠. For moderately large $N$⁠, one will see deviations between the trajectories that are due to the stochastic fluctuations of the ABM dynamics. These fluctuations can be included into the description by an extension of the ODE system to a stochastic differential equation (SDE),^20,21 which allows us to reproduce the fluctuations up to second order. Since the ABM process and the SDE process are of stochastic nature, it is often omitted that this is a pathwise result for medium to large $N$ and fixed finite time interval $[0,T]$⁠, too, as we will also discuss in more detail below. In most articles on the subject, the SDE-based description is not discussed directly, but instead in terms of a mean-field approximation leading to the corresponding Fokker–Planck equations (see, e.g., Refs. 8 and 22).

When there is dynamical behavior on timescales that scale exponentially with $N$⁠, i.e., $T≍exp⁡(ζN)$⁠, case (B) from above, formal approaches via asymptotic expansions using the well-known WKB method²³ or by means of the mathematically more rigorous large deviation principle (LDP)²⁴ allow us to accurately characterize the associated rare events as remote tail events. That is, there are rare events that appear with very small probability $p≍exp⁡(−ζN)$ as large deviations from the expected dynamical behavior of the ABM. As one example, we consider ABM dynamics, which exhibit so-called metastable sets or states to which the process is attracted for long periods of time before the rare event of an exit from the metastable sets and into another one occurs. For social systems, understanding such rare events can be of great interest because they should be avoided (e.g., an outbreak of pandemics or riots) or because, on the contrary, a transition is desirable (e.g., a shift to environmentally friendly technologies or healthy lifestyles). Often, the respective rare events, e.g., switching between metastable sets, cannot be described by the mean-field ODE, which typically exhibits asymptotically stable fixed points instead of metastable sets. Thus, the ODE system often is irrelevant for studying the behavior of the ABM on the longest timescales. As we will see below, the rare events can be understood and characterized by LDP-based analysis, at least in theory. However, for large $N$⁠, the necessary computations often are infeasible in practice. As has been observed and discussed in the literature,²⁵ in some cases, they can be replaced by LDP-based computations for the associated mean-field SDE, at least under certain precautions that will be discussed below in detail.

The connection between (A) and (B) lies in understanding the behavior of the system on long (but not necessarily exponential) timescales connected to metastability. For this case, several approaches utilizing transfer operators have been developed, e.g., so-called Markov state models (MSMs).^26,27 We will show how to adapt and use them for studying ABM dynamics. As the main result, these techniques will allow us to use results from (A) in order to characterize rare events on long timescales for (fixed) medium or large $N$ and to understand how this characterization is related to (B). Similar questions related to metastability have already been discussed in the literature, mainly for (discrete-time) Markov chain approaches,²⁸ and via the Fokker–Planck equation.²⁹ Herein, we will provide new insights illustrated by specific examples. The main point is that MSM building and related techniques like extended dynamic mode decomposition (EDMD)—see Ref. 30 for a recent overview—require many short trajectories of finite length only, allowing, e.g., for parallel computing to determine quantities of interest.

The paper is structured as follows: the agent-based model and its continuous-time description as a transition jump process are introduced in Sec. II. A guiding example reappears in Secs. III–V for illustration. We rescale the transition jump process and consider the large population limit processes given by ODE and SDE systems. Next, long timescales and metastability are discussed and illustrated in Sec. III by means of the transfer operator approach. Using MSMs and the projected transfer operator, we characterize mean first exit times. These reappear in Sec. IV, where we analyze rare events on exponentially long timescales using LDP-based approaches. We address the advantages and limitations arising in both approaches. Finally, open questions and future work are discussed in Sec. V.

The starting point in this paper is a continuous-time ABM without spatial resolution where agents are nodes in an interaction network, similarly to the discrete-time model given in Ref. 31. In analogy to the term interacting particle system, which is used in the context of chemical reaction networks, we also refer to the considered model as an interacting agent system. Each agent can switch according to given transition rules between finitely many states (called types). If the type space is ${1,…,n}$ and we consider $N$ agents, then the state space of the ABM is ${1,…,n}N$⁠, i.e., its size grows combinatorially like $nN$⁠, which poses a problem when $N$ becomes very large.

Alternatively, the ABM can be described via the population state, which counts the number of agents for each of the available types. That is, the size of population state space scales asymptotically like $Nn$ in the worst case. For the basic case of homogeneous agents and a complete interaction network, or if agents randomly interact with each other, the population space description is exact, otherwise it already entails an approximation due to aggregation. In the following, we assume completeness of the interaction network for the sake of simplicity. Then, the ABM transition rules between the types imply transition rules between population states. The associated propensities depend on the total population state.

The temporal evolution of the population state is given by the Markov jump process, called transition jump process or ABM process, and can be characterized by a master equation, in direct analogy to the reaction jump process considered in the context of chemical reaction systems. By applying Gillespie’s algorithm³² for the numerical simulation of the transition jump process, the problems of scheduling described above are naturally avoided. However, these simulations become computationally intensive for settings with numerous agents.

The considered system of interacting agents consists of

• a fixed number $N∈N$ of agents,

• a set ${Si:i=1,…,n}$ of types $Si$ available to the agents, i.e., at any point in time each agent owns one of these types,

• a set ${R1,…,RK}$ of transition rules $Rk$ defining possible changes of the agents’ types, i.e., representing actions of and interactions between agents,

• a set of propensity functions specifying the rates at which transitions randomly occur depending on the population state.

In the agent-based context, transition rules are mostly described from the perspective of a single agent; they may also be referred to as update rules or decisions. To completely specify the dynamics, in this case, also a description of the interaction patterns and update orders (scheduling) is needed. Here, we adopt the formulation used commonly in the chemical context that takes a system level perspective, leading to a simpler and more transparent description: each transition rule $Rk$ is represented by an equation of the form

with the coefficients $alk,blk∈N0$ denoting the numbers of agents of each type involved in the transition as input and output interaction partners (see Example II.1). In order to ensure that the total number of agents is not changed by such a transition, we assume $∑i=1naik=∑i=1nbik$ for each $k$⁠. The associated vector $νk=(ν1k,…,νnk)∈Zn$ is defined by

and describes the net change in the number of agents of each type $Si$ due to transition $Rk$⁠. Note that different transitions can lead to identical net changes, accounting, e.g., for different social mechanisms inducing a specific change. The population state of the system (i.e., the system state) is given by a vector of the form

where $xi$ refers to the number of agents of type $Si$⁠. As the total number of agents is assumed to be constant, the population state space³³ is given by

Each transition $Rk$ induces an instantaneous change in the population state of the form

which may occur at any point in time $t≥0$⁠. The probability for transition $Rk$ to happen in an infinitesimal time step $dt$ is given by $αk(x)dt$⁠, where $αk:X→[0,∞)$ is the propensity function of this transition. We assume the propensity $αk(x)$ to be proportional to the number of combinations of interacting agents in $x$⁠, and, moreover, to scale with the total population size $N$⁠.

One way to describe the temporal evolution of the system is to formulate a continuous-time stochastic process $X(t)|t≥0$⁠,

with $Xi(t)$ denoting the number of agents of type $Si$ at time $t≥0$⁠. The process $X(t)|t≥0$ is a Markov jump process, i.e., it is piece-wise constant, with jumps of the form

occurring at random jump times. The waiting times between the jumps follow exponential distributions determined by the propensity functions. Let $P(x,t):=P(X(t)=x∣X(0)=x0)$ denote the probability to find the system in state $x$ at time $t$ given some initial state $x0$⁠. In order to describe the time-evolution of the system, we consider the Kolmogorov forward equation of the transition jump process $X(t)|t≥0$⁠, which is given by

where we set $αk(x):=0$ and $P(x,t):=0$ for $x∉N0n$ in order to exclude terms in the right-hand side of (2) where the argument $x−νk$ contains negative entries. In the context of chemical reaction networks, (2) is called chemical master equation.³⁷ 

Alternatively, the evolution of the probability $P(x,t)$ can be described by the transfer operator $Pt$ acting on functions of $x$ as

That is, the transfer operator is the evolution operator $Pt=exp⁡(tG)$ generated by the master equation $dP/dt=GP$ with $G$ denoting the operator on the right-hand side of (2). Generally, the master equation (2) cannot be solved analytically. Instead, the process’ distribution is typically estimated by Monte Carlo simulations of the underlying Markov jump process $X(t)|t≥0$⁠, which may be generated using Gillespie’s stochastic simulation algorithm.³² With an increasing population size $N$⁠, the jumps in the population process $X(t)|t≥0$ occur more and more often, while the size of an individual jump relative to the agent numbers $xi$ decreases. In particular, in social systems, it is also obvious that in a larger system, more actions take place and that generally the influence of the single agent diminishes. Stochastic simulations of the jump process using Gillespie’s algorithm become inefficient for large $N$⁠, since they track each of the individual stochastic jump events, and thus time advance in each iteration becomes extremely small. Thus, simulation of the ABM process may become computationally infeasible for very large values of $N$⁠.

For large $N$⁠, an appropriate rescaling of the dynamics may lead to effective approximate dynamics given by stochastic or ordinary differential equations. In what follows, we summarize the corresponding results, which are mainly based on the law of large numbers and the central limit theorem.

Let $XN(t)|t≥0$ denote the transition ABM process as defined in Sec. II B given the total number $N$ of agents. Also, for the propensity functions, we now use the notation $αkN$ in order to indicate their dependence on $N$⁠. We introduce the relative frequencies $c=x/N$ and rewrite the master equation (2) of the interacting agent system in terms of the frequency-based probability distribution

together with the propensities

where $ε:=N−1$ denotes the smallness parameter. Then, with adapted notation, the master equation (the forward Kolmogorov equation) reads

The transfer operator associated with the rescaled process will be denoted by $P~Nt$⁠, i.e., $P~Nt=exp⁡(tG~N/ε)$ where $G~N/ε$ denotes the generator of the rescaled master equation (3).

It is well-known that in the limit $N→∞$ and assuming convergence of the scaled propensity functions $α~kε⟶ε→0α~k$⁠, the rescaled jump process $XN(t)/N|t≥0$ converges to the frequency process $C(t)|t≥0$ given by the ODE,

with initial state $C(0)=limN→∞XN(0)/N$⁠.³⁸ The order of convergence is given by $N−1/2$⁠,^20,38 i.e., we have the pathwise approximation

for every finite $T$ and an a.s. finite constant $ζ$⁠, both independent of $N$⁠, where $ζ$ will, in general, depend on $T$⁠. Here, $∥⋅∥$ refers to the Euclidean norm (or any other norm) on $Rn$⁠.

A higher order approximation can be obtained by considering the stochastic Langevin process $CL(t)|t≥0$ given by the SDE,

which, in the context of chemical reaction kinetics, is known as the chemical Langevin equation.³⁹ Here, $Bk(t)$⁠, $k=1,…,K$⁠, are independent Brownian motion processes. For the pathwise approximation error, one gets from^20,21 that

for finite $T$ and an a.s. finite constant $ζ$⁠, both independent of $N$⁠, again as above. Also, the estimate is a pathwise result again. It is valid if in (5), one uses a realization of the Brownian motion that is adapted to the realization of the rescaled Poisson process used for the jump process, for details visit.²⁰ Note that these approximation results only hold for finite time intervals that are independent of $N$⁠.

The Fokker–Planck equation associated with the SDE limit process (5) has the form

with

(see Ref. 8 for the discrete-time case). The forward transfer operator $TNt$ associated with the SDE limit process again is the evolution operator

generated by the Fokker–Planck equation $∂/∂tρL=GNρL$⁠, that is, by the operator $GN$ on the right-hand side of Eq. (7). Just as the master equation of the Markov jump process, also the Fokker–Planck equation (7) in general cannot be solved analytically. However, trajectories of the Langevin process may be generated using the standard Euler–Maruyama scheme. Here, possible stepsizes are asymptotically independent of the population size $N$⁠, which makes the model more efficient than the jump process model if the number of agents is large.

Under some technical conditions (e.g., regarding growth of $b$ and $Σ$ at infinity), a unique invariant measure (stationary density) exists, however, perhaps only on the accessible part of the state space. The SDE dynamics is reversible (i.e., satisfies detailed balance with respect to this measure) under some technical conditions (e.g., growth at infinity), if and only if there exists a smooth function $V$ such that⁴³ 

with the more explicit form

If a solution $V$ with sufficient growth at infinity exists, the invariant measure $μ^$ of the SDE dynamics (more exactly the density associated with it) is given by

where $Z$ is a normalization factor.

In order to motivate the investigations in Sec. III, we now calculate the invariant measure for the guiding example depending on $N$⁠, observing that its properties cannot fully be explained by large population limits.

For appropriately chosen rate constants $γij$ and $γij′$⁠, the ABM process of the guiding example given in Sec. II A can exhibit metastable behavior, where one type predominates the population for a long period of time before another type will take over eventually and then dominates for another long period of time. In many (but not all) cases for moderate to large $N$⁠, such metastable behavior can be reproduced by the SDE limit process, while the ODE limit process often fails to reproduce it. As finite timescales (i.e., independent of $N$⁠) may not be sufficient to understand this, we first introduce the transfer operator approach to metastability and then show how it can be applied to the ABM and limit SDE process.

In the following, we will consider the transfer operators $P~Nt$ of the transition jump (ABM) process and $TNt$ of the corresponding limit SDE, where the subindex $N$ indicates for which number of agents the respective dynamics is considered.

For the next step, we assume that both the ABM dynamics and the SDE limit dynamics are (geometrically) ergodic on the respective accessible state space, and there are respective invariant measures $μ$⁠, respectively $μ^$⁠. Moreover, we assume that both kinds of dynamics are reversible. This assumption makes most of the following constructions much shorter and easier to comprehend because, for reversible dynamics, the associated transfer operators are self-adjoint in the Hilbert space $Lπ2$ (with $π=μ$ or $π=μ^$⁠, respectively) and thus possess real-valued spectral decompositions. The assumption of reversibility is not necessary; without it, one has to replace real-valued spectral decompositions by the more involved complex ones, alternatives like the Schur decomposition⁴⁴ or the respective singular value decompositions.

Let $λlt$ be the eigenvalues of $P~Nt$⁠, sorted by decreasing absolute real value, and $ψl$ the corresponding eigenfunctions, where $l=1,2…$⁠. Under our assumptions, it holds that $λ1=1$ is independent of $t$⁠, isolated and the sole eigenvalue with absolute value $1$⁠. Furthermore, $ψ1=μ$⁠. The subsequent eigenvalues satisfy $λlt=exp⁡(−νlt)$ for some $νl>0$⁠, $l=2,3,…$⁠, and thus decrease monotonously to zero both for increasing index and time. That is,

The associated eigenfunctions $ψ2,ψ3,…$ can be interpreted as sub-processes of decreasing longevity in the following sense: Given a function $u$ and lag time $τ>0$ with $u=∑l=0∞βlψl$⁠, $βl∈R$⁠, then

since there exists an index $d∈N$ such that $λlt≈0$ for all $t≥τ$ and all $l>d$⁠. Hence, the major part of the information about the long-term density propagation of the ABM dynamics is encoded in the $d$ dominant eigenpairs. The timescales associated with the dominant eigenvalues $λlt$⁠, $l=2,…,d$ are called dominant timescales and given by

where $Tl$ is independent of $t$⁠, since the transfer operator exhibits an infinitesimal generator, and, therefore,

For the SDE transfer operator $TNt$⁠, we denote the eigenvalues by $λ^lt$⁠, and the eigenfunctions by $ψ^l$ with $ψ^1=μ^$⁠. The same relation concerning the dominant timescales $T^l$ holds.

The ABM transfer operator $P~Nt$ acts on functions $u∈Lπ2(X/N)$ on the rescaled population space $X/N$⁠, while $TNt$ acts on functions $u∈Lπ2([0,1]n)$⁠. In general, the size of these function spaces prohibit efficient computation of the dominant eigenvalues. Therefore, one considers projected transfer operators instead.

We make the following considerations for the transfer operator $P~Nt$ of the ABM process; however, analogous statements hold for the transfer operator $TNt$ of the SDE process. The projected transfer operator with respect to a subspace $L⊂Lπ2(X/N)$ is the operator $QP~NtQ$⁠, where $Q:Lπ2(X/N)→L$ denotes the orthogonal projection of the original function space of $P~Nt$ to $L$ with respect to the appropriate scalar product $⟨⋅,⋅⟩π$⁠. If $L$ is finite-dimensional and spanned by the functions $ϕi$⁠, $i=1,…,m$⁠, then

The most prominent case is that of subspaces spanned by indicator functions $ϕi=1i$⁠, where $1i$ is one on a set $Bi$ and zero outside and the $Bi$⁠, $i=1,…,m$⁠, form a complete, non-overlapping partition of the respective state space. In this case, called the box discretization of the transfer operator, the projected transfer operator $QP~NtQ$ has a matrix representation denoted $Pt=(Pij)i,j=1,…,m$⁠, with entries

where the subindex $π$ means that initial values are distributed due to the respective invariant measure $π=μ$⁠, or $π=μ^$⁠, respectively, and $c$ the respective ABM or SDE process. The matrix $Pt$ representation of the projected transfer operator is stochastic. Its entries can be approximated by starting $n0$ trajectories of the respective process in each box. Let $cik(0)$⁠, $k=1,…,n0$ be initial values distributed due to $π|Bi$ and $cik(t)$ the final point of a $t$-long realization of the process started in $cik(0)$⁠, then

There are a variety of results on how closely the dominant eigenvalues of the projected transfer operator $QP~NtQ$ approximate the eigenvalues of the original transfer operator $P~Nt$⁠: for example, in Ref. 45, the error is characterized via the projection error induced by $Q$⁠, while in Ref. 26 (Thm. 4.16), upper and lower bounds to the error, specific for every complete partition of the respective state space into sets, are given. Moreover, in Ref. 46, the long-term error $∥P~Nst−(QP~NtQ)s∥$ with $s≫1$ is analyzed.

After these preparations, we can point out the key idea of the transfer operator approach to metastability:²⁶ metastable subsets can be detected via the dominant eigenvalues of the respective transfer operator close to its maximal eigenvalue $λ=1$⁠. Moreover, they can be identified by exploiting the corresponding eigenfunctions $ψl$⁠, $l=1,…,d$⁠. In doing so, the number of metastable subsets is equal to the number of eigenvalues close to 1, including 1 and counting multiplicity. Weighted with the invariant measure $π$⁠, the dominant eigenfunctions $ψl~=ψl/π$ exhibit different sign combinations for each set (see Ref. 26). Alternatively, they allow us to decompose state space into metastable sets by their zeros.⁴⁷ Then, the asymptotic exit rates from these metastable sets are approximately given by the inverse of the dominant timescales.⁴⁷ 

The dominant eigenfunctions also play a major role in the theoretical justification of building so-called Markov state models (MSM).^26,27 An MSM is a very low-dimensional reduced model of the full transfer operator $P~Nt$ with almost the same dominant eigenvalues. For $d$ dominant eigenvalues, the Markov states of an MSM are formed by non-negative ansatz functions $ϕk$⁠, $k=1,…,d$ that satisfy $P~Ntϕk≈ϕk$ as closely as possible and thus can be interpreted as macro-states that are almost invariant under the dynamics. More precisely, the $ϕk$ are selected such that they form a partition of unity, i.e., $∑kϕk=1$⁠, and the projected transfer operator $QP~NtQ$ (represented by a $d×d$ matrix) satisfies the condition that its eigenvalues are as close as possible to the dominant eigenvalues of $P~Nt$⁠. There are several algorithms for computing these almost-invariant ansatz functions $ϕk$ based on the dominant eigenfunctions of $P~Nt$⁠, e.g., set-oriented algorithms like milestoning,⁴⁸ or algorithms approximating a basis of the dominant eigenspace like PCCA+⁴⁹ (see Example III.1 for more details).

All the previous considerations can be carried out for both the ABM process and the SDE process (replacing $P~Nt$ by $TNt$⁠), independently of $N$⁠. Moreover, for large enough $N$⁠, the respective results of this MSM analysis are very close for the ABM and SDE cases and converge for $N→∞$⁠, which we will see in the following two examples.

In order to get additional information, we are interested in the expected mean first exit time $ηN(i)$ of the two processes for starting in some box $Bi$⁠, $i∈I$ with $I={i=1,…,m:Bi∩A=∅}$ and going into some set $A$⁠. Note that $A$ and $Bi$ do not need to be metastable sets. The vector $ηN=(ηN(i))i∈I$ of these mean first exit times can be computed by solving

where $PIτ$ denotes the submatrix of the $m×m$ discretization matrix $Pτ$ of the projected transfer operator $QP~NτQ$ (respectively, $QTNτQ$⁠) for indices $i∈I$⁠, and $1I$ the vector of length $|I|$ with all entries having value one. Again, Eq. (11) results from the Galerkin discretization of the respective equation of the full transfer operator.

Given the mean first exit time $ηN$⁠, we can also compute the vector $ΦN=(ΦN(i))i∈I$ of associated rates via

We will calculate and compare the mean first exit times in the following example.

The closeness between the discretization matrices $PSDEt=QTNtQ$ and $PABMt=QP~NtQ$ of the transfer operators of the SDE and ABM processes for large $N$ is due to the pathwise closeness (6) of the SDE and ABM process, which implies

entrywise^20,26 for all times $t$ in $[0,τ]$ for not too large $τ$⁠. If the ansatz space used for the discretization is of dimension $m$⁠, then this means that there is a constant $ζ0<∞$ such that

as an asymptotic result, i.e., there is an $N0∈N$ such that the statement holds for all $N>N0$⁠. This result requires only trajectories of length $τ$⁠, with the additional computational advantage that they can be computed in parallel for all boxes. This implies that for large enough $N$⁠, we may characterize the long-term behavior of the ABM process via the discretized SDE transfer operator.

However, the discretizations $PSDEt$ and $PABMt$ have to be computed by means of trajectories via (10) (denoted by $PSDE,n0t$⁠), which introduces an additional sampling error of order $1/n0$ in the number $n0$ of trajectories. Moreover, the discretizations differ from the true transfer operators by the discretization error $errdiscr$ resulting from using a finite-dimensional ansatz space (see Ref. 45,46) such that in total

where $P~Nt$ denotes the true ABM transfer operator. Note that this estimate does not require any assumption of a gap in the spectrum of $P~Nt$ (as often is assumed inappropriately) (see Ref. 50 and, in particular, the discussion in Ref. 46).

Therefore, even for large $N$⁠, metastable structures of the underlying ABM process can be identified by studying $PSDEt$ only up to the resolution given by the underlying discretization, which may be limited according to computational resources. In particular, if the metastable structures became finer and finer for increasing $N$⁠, one would have to increase resolution by increasing the dimension $m$ of the ansatz space with $N$ in order to yield sufficient accuracy which, in turn, would mean that we might not be able to reduce the first error term in Eq. (13) sufficiently.

The approximation results of Sec. II C concern the pathwise behavior of the stochastic ABM process on finite time intervals (that do not scale with $N$⁠). In Sec. III, we went beyond finite timescales but saw that some rare events associated with metastability happen on timescales that grow exponentially with $N$⁠. In this case, i.e., when $T≍exp(ζN)$⁠, the transfer operator approach reaches its limitations due to the increase of dimension $m$ of the ansatz space (see Sec. III B). We now investigate these unlikely tail events of the dynamics and their approximation by the SDE process via the LDP.

As for the guiding Example II.1, we will see that (under specific conditions outlined below) there is an approximate agreement of the tail probabilities when comparing the ABM process to the corresponding SDE process. The literature contains results that show that this is not always the case (see, e.g., Ref. 23). We will shed some light on the reasons behind success and failure too. The analysis is done in terms of the rescaled ABM process $Cε(t):=1NXN(t)$ for $ε=1/N$⁠, with the corresponding master equation given in (3). Again, $ε=1/N$ denotes the smallness parameter for the sake of keeping the notation that is most prominent in large deviation theory, for which we will now give a brief introduction.

Let $Cε(⋅)=Cε(t)|t∈[0,T]$ be an (arbitrary, random) path from an appropriate path space $C$ [e.g., $C=H1([0,T],Rn)$] starting in $Cε(0)=c0$⁠, and let $P$ denote the probability distribution generated by the master equation (3), respectively, by the Markov jump process associated with it, or by the SDE limit Eq. (5). Furthermore, let $c(⋅)=c(t)|t∈[0,T]$⁠, $c(0)=c0$⁠, denote a specific path from the selected path space $C$ [for example, the solution trajectory of the ODE limit system (4)]. Then, we call $I:C→[0,∞]$ the large deviation rate function associated with the law $P$ if

(see Ref. 54). This is often written in the following, more intuitive form:

where $≈$ denotes $δ$-closeness of the paths $Cε(⋅)$ and $c(⋅)$⁠, and $≍$ asymptotic equality or exponential equivalence, that is,

where the notation $xε=x+o(1)$ means that $(xε−x)/ε→0$ for $ε→0$⁠. The statement says that the probability to find a solution path $Cε(⋅)$ that is close to the specific path $c(⋅)$ is exponentially small in $ε$ with an asymptotic rate given by the rate function $I$⁠. In an alternative way of presenting this, $I$ is the path space measure introduced by the dynamics for small $ε$⁠: Following Ref. 55, we can write the probability $pε(t0,c0,t1,c1)$ to go from $c0$ at time $t0$ to $c1$ at time $t1$ using the path integral formalism:

where $Dc(⋅)$ denotes integration over all paths $c(⋅)=c(t)|t∈[t0,t1]$ with $c(t0)=c0$ and $c(t1)=c1$⁠. That is, for small $ε$⁠, the exponential factor $exp⁡(−1εI(c(⋅)))$ works as a probability density in the space of paths.

In many cases, expressions for the rate function can be found by constructing its point-wise form $I:Rn×[0,T]→[0,∞]$⁠, such that the intuitive Eq. (14) becomes

meaning $εlogP(Cε(t)≈c)=−I(c,t)+o(1)$ for all $t$⁠. While $I$ characterizes the asymptotic behavior of the probability distribution $P$ induced by the dynamics in the state space $Rn$⁠, $I$ characterizes the associated probability distribution in the path space $C$⁠. In what follows, we will see how $I$ can be computed for the ABM process (3) and for the SDE (5). We will present just the essential results; there are rigorous construction techniques for the rate functions like the Feng–Kurtz⁵⁴ or the Gärtner–Ellis⁵⁶ methods (see Ref. 24). Furthermore, there are several other asymptotic techniques for studying exponentially small probabilities for master and Fokker–Planck equations in the sense of Eq. (15), e.g., WKB theory or eikonal approximations (cf. Ref. 23,25).