Dynamic maximum entropy provides accurate approximation of structured population dynamics
DDYNAMIC MAXIMUM ENTROPY PROVIDES ACCURATE APPROXIMATION OFSTRUCTURED POPULATION DYNAMICS
Katar´ına Bod’ov´a a ∗ , Enik˝o Sz´ep b , Nicholas H. Barton ba Faculty of Mathematics, Physics and Informatics, Comenius University,Mlynsk´a Dolina, 84248 Bratislava, Slovakia, b Institute of Science and Technology Austria (IST Austria),Am Campus 1, Klosterneuburg A-3400, Austria ∗ [email protected] Abstract.
Realistic models of biological processes typically involve interacting components on multiplescales, driven by changing environment and inherent stochasticity. Such models are often analytically andnumerically intractable. We revisit a dynamic maximum entropy method that combines a static maximumentropy and a quasi-stationary approximation. This allows us to reduce stochastic non-equilibrium dynamicsexpressed by the Fokker-Planck equation to a simpler low-dimensional deterministic dynamics, without theneed to track microscopic details. Although the method has been previously applied to a few (rathercomplicated) applications in population genetics, our main goal here is to explain and to better understandhow the method works. We demonstrate the usefulness of the method for two widely studied stochasticproblems, highlighting its accuracy in capturing important macroscopic quantities even in rapidly changingnon-stationary conditions. For the Ornstein-Uhlenbeck process, the method recovers the exact dynamicswhilst for a stochastic island model with migration from other habitats, the approximation retains highmacroscopic accuracy under a wide range of scenarios for a dynamic environment. Introduction
Conceptual understanding of realistic problems in applied sciences is often hindered by the curse ofcomplexity, with quantities of interest coupling to finer features. Due to their multiscale character evensimple questions lead to exploration of the full complexity of the system. But how can we ever understandthe processes around us if incremental learning is impossible?Statistical mechanics provides a clever way to understand complex multiscale problems by linking processeson different scales through the parsimony principle – the method of maximum entropy (ME), introduced by[1]. ME has the form of a variational problem where an entropy of the microscopic distribution is maximized,while enforcing macroscopic constraints, e.g. average energy of gas particles [1], see figure 1A. The methodgained popularity in applied sciences in recent decades primarily as a tool for inference from empirical data,for instance in bird flocking [2], neuronal firing [3], or protein variability [4].However, realistic biological questions often do not adhere to the assumption of stationarity. Adaptationof populations to spatial and temporal ecological gradients is an example of a complex non-equilibriumprocesses in which ecological and evolutionary processes interact [5–9]. How can a simple concept of ME beapplied to such systems? The most straightforward approach is to use ME on dynamic trajectories, forcingconstraints on the dynamical features (figure 1B). This approach, called maximum caliber (MC), introducedby [10] is suitable for inference of models and their parameters from dynamic data. Comprehensive reviewsof MC [11–13] provide multiple examples where the method has been successfully applied to non-stationarybiological processes.While such an inverse approach is useful for understanding temporal data, in many cases data are notavailable but the dynamics, although often extremely complicated, are known to be accurately describedby the Fokker-Planck equation (FPE). The aim of our work is to demonstrate usefulness of a theoreticaldimensionality-reduction technique, which approximates possibly many-dimensional stochastic dynamics toa low-dimensional deterministic dynamics of a few key observable quantities (figure 1C). The approximationcombines static ME ansatz in the FPE equation with a quasi-stationary assumption. We refer to this method a r X i v : . [ q - b i o . P E ] F e b A)Maximum entropy (static approximation)Constraints, i.e. time time (C)Dynamic maximum entropy
Combination of MaxEnt and quasi-stationarityDynamics of effective forces: time (B)Maximum caliber method (MaxEnt distribution over trajectories)Constraints, e.g.
Figure 1.
Variational methods ME and MC compared to DME. (A) ME looks at a snapshot x of a process at a particular time and provides an approximation ¯ u ME ( x ) of the microscopicdistribution, given knowledge of a few key macroscopic observables. (B) MC is analogousto ME, however, each data point represents a trajectory x ( t ). MC connects the microscopicdistribution over possible trajectories with macroscopic constraints and approximates it by¯ u MC ( x ( t )). (C) DME is a quasi-stationary approximation of the stochastic dynamics, givenby the FPE, which reduces the full problem to a low-dimensional dynamics. This reductionis a consequence of a ME ansatz; the approximation at each time ¯ u DME ( α ( t )) solves theME problem (stationary form in the FPE), where the dynamics of the effective forces α aresystematically derived from the FPE.as the dynamic ME (DME). DME has been used before to solve problems in quantitative genetics [14–16]and independently in cosmology [17–20].The most surprising feature of DME is its accuracy. The method, derived from the assumption of quasi-stationarity, remains extremely accurate even in far-from-equilibrium regimes [14–16, 21]. Nevertheless,explicit estimates of the method’s accuracy opens a challenging mathematical problem. While in reactionkinetics, where the quasi-stationary assumption is often used to reduce dynamics of reactants by assumingthat concentration of certain chemicals does not change at the considered timescale, its validity has beenshown rigorously only for a few basic systems using singular perturbation [22–24]. In comparison, theaccuracy of DME still remains a mystery. We will address this puzzle by focusing on two processes: theOrnstein-Uhlenbeck (OU) process, for which DME gives an exact solution, and a logistic model of populationgrowth in a continent-island model, which is one of the simplest population models that includes geographicpopulation structure. These processes are explored in the far-from equilibrium regime.2. Dynamic maximum entropy
Here we present a dynamic maximum entropy (DME) method to approximate stochastic dynamics in termsof a FPE [14–21]. The method is based on a combination of ME in statistical physics [1], which solves thestationary problem exactly, with a quasi-stationary assumption, as typically used in chemical kinetics [22]to reduce the number of equations. The method applies to stochastic dynamics with an explicit stationarydistribution, even though its application is not limited to such problems (as shown in [21] a solution ansatz,which is not based on the stationary form can sometimes lead to more accurate approximation). DMEwas introduced in population genetics to understand how quantitative traits change in time in the presenceof various evolutionary mechanisms without resolving details about the dynamics of the underlying genefrequencies. Independent use of the method in statistical physics focused on exact and approximate solutionsof the nonlinear FPE arising in cosmology. We provide a comprehensive summary of the method based on[21] in this section. ssume stochastic dynamics in the Langevin form(1) d x k ( t ) = g ( x k )2 ∂∂x k (cid:34) d (cid:88) i =1 α i A i ( x ) (cid:35) d t + g ( x k )d ξ k ( t ) . for x = ( x , . . . , x N ), x k ∈ Ω X , t >
0, and x (0) = x . The potential in the first term is a linear combinationof forces α i , acting on functions A i ( x ), which may introduce coupling between equations. Function g ( x k )represents amplitude of stochastic fluctuations and ξ k ( t ) are independent Wiener processes. Previous studies(e.g. [14–16, 25]) focused on examples in population genetics where x k corresponds to the frequency of acertain gene, affecting some quantitative trait. This frequency depends on evolutionary processes, e.g.selection, mutation, and inherent stochastic fluctuations, described by the forces α i . At this point weconsider constant forces α i ∈ R , which means that the distribution u ( t, x ) follows dynamics described by theFPE(2) ∂u ( t, x ) ∂t = − N (cid:88) k =1 d (cid:88) i =1 α i ∂∂x k (cid:20) g ( x k )2 ∂A i ( x ) ∂x k u ( t, x ) (cid:21) + 12 N (cid:88) k =1 ∂ ∂x k (cid:2) g ( x k ) u ( t, x ) (cid:3) , which can be also expressed in the flux form ∂ t u ( t, x ) = − (cid:80) Nk =1 ∂ x k J k [ t, x ]. This FPE is complementedwith no-flux boundary conditions J k [ t, x ] = 0 at x k ∈ ∂ Ω X and the initial condition u (0 , x ) = u ( x ).2.1. Stationary solution and ME.
At large times the dynamics approach a stationary distribution,parametrized by the vector of forces α = ( α , . . . , α d ) ∈ R d (3) ¯ u α ( x ) = 1 Z α (cid:32) N (cid:89) k =1 g ( x k ) (cid:33) exp (cid:34) d (cid:88) i =1 α i A i ( x ) (cid:35) with the normalization coefficient (i.e., the partition function)(4) Z α = (cid:90) Ω NX (cid:32) N (cid:89) k =1 g ( x k ) (cid:33) exp (cid:34) d (cid:88) i =1 α i A i ( x ) (cid:35) d x , where A = ( A , . . . , A d ) is a vector function of the state variables x which the forces α i acts on. Using theterminology in statistical physics we refer to functions A i as observables, as their expectations in problemsin physics provide a macroscopic description of the system in terms of its natural observable quantities (i.e.,average energy of a gas particle, as formulated by [1]). We extend our scope from looking at a problemwith constant forces α to time-dependent forces α ( t ) to account for realistic scenarios where the dynamics,initially settled to a stationary solution, are pushed out-of equilibrium by changes in the forces α . Thedynamics of the expectation (cid:104) A j (cid:105) follows from the FPE. Using notation B ji = (cid:68)(cid:80) Nk =1 g ( x k )2 ∂A j ∂x k ∂A i ∂x k (cid:69) , V j = (cid:68)(cid:80) Nk =1 g ( x k ) ∂ A j ∂x k (cid:69) we obtain ∂∂t (cid:104) A j (cid:105) = d (cid:88) i =1 B ji α i + 12 V j . (5)This forms a system of ordinary differential equations for (cid:104) A (cid:105) , which is generally not closed due to nonlinearityof the functions A j ( x ). Next we define a logarithmic relative entropy(6) H [ u | ¯ u α ] := (cid:90) Ω NX u ln u ¯ u α d x , where u ( t, x ) is a solution of the FPE at time t ≥ x k ∈ Ω X . For any t ≥ α = α ∗ , which can be obtained by solving a set of first-order conditions(7) 0 = ddα i H [ u | ¯ u α ] = − (cid:90) Ω NX u ¯ u α ddα i ¯ u α d x = (cid:104) A i (cid:105) ¯ u α − (cid:104) A i (cid:105) u . The above intriguing relationship states that if for a given time t there exists a maximum of the relativeentropy (6) with respect to all α i reached for some choice of parameters α ∗ , then the expectation of A i through the distribution u ( t, x ) equals the expectation through the stationary distribution ¯ u α ∗ ( x ) at this ime. This simple fact, shown previously in [21] and in a slightly different form by [14–16] suggests thatinstead of the full representation of the problem using FPE one could trace only the d -dimensional dynamicsof α ∗ , that parametrize the approximate solution by the form (3). Furthermore, the two representationsagree in terms of the expectations (cid:104) A i (cid:105) at a given time. Note that the solvability of the equation (7) (for α ∗ )is nontrivial [21] and requires further attention. Nevertheless, concavity of the relative entropy is implied bythe following relationship(8) d dα i dα j H [ u | ¯ u α ] = (cid:90) Ω NX (cid:20) ddα j ln ¯ u α (cid:21) (cid:20) ddα i ln ¯ u α (cid:21) u d x − (cid:90) Ω NX (cid:20) u α d dα i dα j ¯ u α (cid:21) u d x = Cov( A i , A j ) ¯ u α analogous to similar expressions in [14–16, 21] stating that the Hessian of the relative entropy is positivesemidefinite.2.2. Dynamical approximation.
We have established a relationship between the solution of the full sto-chastic dynamics (2) and a stationary form u α ∗ parametrized by suitable effective forces α ∗ following ME.However, as we demonstrated in figure 1A, ME is applicable only to static problems. When the system isout-of-equilibrium, we need to establish a dynamic relationship between the values α ∗ ( t ) and α ∗ ( t ) for t (cid:54) = t by using the information captured by the FPE.To derive the DME approximation of (2) we use an ansatz u ( t, x ) = ¯ u α ( t ) ( x )+ R ( t, x ) for some continuous α ( t ) where R ( t, x ) is the time-dependent residual. The dynamics of the expectations (5) becomes(9) ∂∂t (cid:104) A (cid:105) α = B α α ( t ) + 12 V α + (cid:20) B R α ( t ) + 12 V R − ∂∂t (cid:104) A (cid:105) R (cid:21) . where (cid:104)·(cid:105) u represents expectation through distribution u and (cid:104)·(cid:105) α := (cid:104)·(cid:105) ¯ u α ( t ) . Now we make two key as-sumptions. First, we assume that the residual terms in the bracket of (11) are small and we neglect them.In addition, we also impose a quasi-stationarity approximation, assuming that α ∗ are chosen to satisfy theequilibrium relationship(10) B α ∗ α ∗ ( t ) + 12 V α ∗ = 0for all t >
0. Both steps are easy to justify if the forces α ( t ) are slowly changing (in the adiabatic regime) andthe solution of the FPE is thus close to an equilibrium form (3). However, its validity when out-of-equilibriumis not clear. We use (10) to replace V α by V α ∗ and B α by B α ∗ in (9) to approximate(11) ∂∂t (cid:104) A (cid:105) α ∗ ≈ B α ∗ ( α ( t ) − α ∗ ) . Finally, to obtain a closed dynamical system for α ∗ we use the chain rule, noting that ∂ (cid:104) A (cid:105) α ∗ ∂t = ∂ (cid:104) A (cid:105) α ∗ ∂ α ∗ ∂ α ∗ ∂t .Combined equations (7) and (8) imply that differentiation of an expectation (cid:104) A i (cid:105) α ∗ with respect to α j givesus a covariance C α ∗ := Cov( A i , A j ) α ∗ . Therefore(12) ∂ α ∗ ∂t = C − α ∗ B α ∗ ( α ( t ) − α ∗ ) , α (0) = α together with a parametric form (3) represents DME approximation of dynamics (2). Solution of (12) can beplugged into the stationary parametric form (3), which allows us to not only study the accuracy of the keymoments used in DME, but also to compute any statistical feature of the approximate solution and compareit with the exact solution.Note that the equation (12) can be solved for any prescribed continuous function α ( t ) and in the mostextreme case, the forces α ( t ) can contain step changes, which clearly violate the quasi-stationarity assump-tion. However, all previously studied applications of DME approach showed that the approximation capturesextremely well the expectations of the key functions even when the forces change rapidly. This is one of themost remarkable and unexpected features of the DME approach which will be studied here.Realistic situations (e.g. selection acting on quantitative traits that depend on many alleles) typicallyinvolve high-dimensional stochastic dynamics with nonlinearities and coupling terms [14–16]. However, forsimplicity we consider examples leading to a simpler one-dimensional form with a remark that more complexmodels can be analyzed using this approach as long as the stationary distribution of the FPE is explicit. owever, even when x is a scalar A and α are vectors in all problems studied here: these vectors summarisethe infinite-dimensional distribution of x .3. Ornstein-Uhlenbeck process
The model.
Here we outline a simple example of stochastic dynamics where DME reproduces the exactdynamics. We use an OU process, which describes dynamics with linear relaxation to an equilibrium in thepresence of constant Gaussian noise (examples include a particle under friction, animal motion, financialtime series, etc.). The OU process has three parameters: µ – long-time average of the state variable x , β –persistence length, and σ – magnitude of noise. It has the form(13) dx = β ( µ − x ) dt + σdξ ( t ) . The stationary distribution of the equation (13) can be obtained from the FPE, which describes time-evolution of the probability distribution of x , denoted by u ( t, x )(14) ∂∂t u = − ∂∂x [ β ( µ − x ) u ] + σ ∂ ∂x u , by setting the left-hand side equal to 0(15) ¯ u α ( x ) = 1 Z exp (cid:20) µβσ x − βσ x (cid:21) = 1 Z exp (cid:20) σ α · A (cid:21) , where Z = 1 / (cid:112) β/ (2 πσ ) is the normalization factor, α = (2 µβ, − β ) and A = ( x, x ). If we set theinitial condition of (13) as the stationary distribution corresponding to the parameters ( β , µ , σ ), i.e., asa Gaussian x ∼ N ( µ , σ / β ), then the solution of (14) is a Gaussian N ( m ( t ) , v ( t )) for every t ≥
0. Thefunctions m ( t ) and v ( t ) represent a time-dependent mean and variance, which solve the system of ODEs˙ m = β ( µ − m ) , ˙ v = σ − βv , (16)with initial conditions m (0) = µ and v (0) = σ β . The explicit solution m ( t ) = µ + ( µ − µ ) e − βt , (17) v ( t ) = σ β e − βt + σ β (1 − e − βt )(18)satisfies m ( t ) → µ and v ( t ) → σ / β as t → ∞ .The stationary solution ¯ u α ( x ) in (15) solves a variational ME problem with the relative entropy definedin (6) using forces α and observables A . Even though the OU process has three natural parameters ( β, µ, σ )ME implies that the stationary solution is a 2-parameter family of functions of form (15). Equivalently, wemay restate the initial distribution by taking a fixed initial value σ = σ and picking β such that the initialvariance of the Gaussian σ / β = v has the desired initial value. This allows us to keep the volatility of theprocess fixed in the DME approach. Following the DME derivation steps (details in Appendix A) we obtaina 2-dimensional dynamical system for the effective forces α ∗ shown in (A.6)-(A.7). This coupled dynamicalsystem can be transformed into decoupled dynamics of µ ∗ and β ∗ of the formd µ ∗ d t = β ( µ − µ ∗ ) , (19) d β ∗ d t = 2 β ∗ ( β − β ∗ ) . (20)First, note that the dynamical system (19)-(20) is independent of the noise magnitude σ , which is consideredconstant. It is explicitly solvable since it is decoupled and while the first equation is linear, the second one islogistic. The DME solution is consistent with (17)-(18) given σ = σ , m ( t ) = µ ∗ ( t ) and v ( t ) = σ / (2 β ∗ ( t ))are functions of effective forces α ∗ . This is due to the linearity of the OU process, which yields a closeddynamics of the first two moments and thus preserves a Gaussian form of the solution at each time, providedwe started with a Gaussian initial condition. Note that the OU process is not the only stochastic processwhere DME provides an exact solution. As [18] showed, the nonlinear extension of the OU process can besolved exactly using a ME ansatz. .2. Numerical example.
Figure 2A shows a numerical simulation of the OU process for three choicesof initial Gaussian distribution (centered at x = 0 . , . , .
2) for 3000 trajectories for each case (allparameters summarized in the figure legend). Initially, the system is in a stationary state corresponding toparameters µ β σ (Gaussian form with parameters x , σ / β ). However, at t = 0 the parameters ofthe OU process rapidly change, pushing the system out-of-equilibrium. As a response, the distribution ofsample trajectories follows a Gaussian form at each time, eventually converging to N ( µ, σ / β ). In figure 2Bwe plot the 2-dimensional DME dynamics of effective forces β ∗ , µ ∗ , which is exact. Each trajectory (vectorfield of the dynamical system is plotted as well) represents the complete solution of the FPE for a givenparameter choice (at each time it is a Gaussian). The microscopic distribution at four different times in panelC shows an agreement between stochastic simulations (histograms) and microscopic distributions obtainedfrom the DME approach (solid curves). In general, the goal of DME is to approximate the dynamics on themacroscale, thus, we do not expect DME to capture also the microscopic distribution. However, for the OUprocess DME is exact and thus the method recovers both the macroscale and the microscale properties ofthe process without loss of precision. Figure 2. (A) Numerical simulations of OU process with parameters β = 0 . µ = 1, σ = 0 .
1. We used three random initial conditions from a distribution N ( x , σ ) with µ = x = 1 . , . , . β = 0 . , . , .
45 and σ = σ . (B) Effective forces ( µ ∗ , β ∗ ) followingdynamics (19)-(20) corresponding to the same set of initial conditions as in panel A. (C)Histograms of x ( t ) at times t = 0 . , , , β = 0 . µ = 0 . σ = 0 . Dynamics of a single island with immigration
Here we use the DME method to approximate stochastic population dynamics. We will study a simple,yet nonlinear model – stochastic logistic population growth in a single island with immigration from otherhabitats. .1. The model.
In the case of unlimited resources and absence of predation, populations would growindefinitely. However, in natural populations this is not the case: various factors impose bounds on thisexponential growth. The logistic growth model describes population size regulation in the absence of de-mographic stochasticity. At low population sizes, when resources are abundant and competition is low, thepopulation grows with its intrinsic growth rate, r . However, the total growth rate of the population decreaseslinearly with increasing population size. In particular, the growth rate is zero when the population is atcarrying capacity, K , reflecting the situation when each individual replaces itself in each generation. Thecarrying capacity represents the maximal sustainable population size. We follow the population in a singleisland, with a migration from other habitats at rate m . In the presence of demographic stochasticity thedescribed population dynamics can be formulated using a stochastic differential equation(21) d n = [ n ( r − λn ) + m ] d t + √ γn d ξ , where n ( t ) represents the population size at time t , λ = r/K is the density regulation and γ describes thevariance in population size. For the sake of simplicity we fix γ = 1 corresponding to a Poisson(1) numberof offspring for each individual (with total variance n ). Extinction in this stochastic dynamics for m = 0 isunavoidable from a mathematical point of view (as the process is a critical branching process) but can beprevented by migration when m > u ( t, n ) (which is of the same form as (2)):(22) ∂u∂t = − ∂∂n [( n ( r − λn ) + m ) u ] + 12 ∂ ∂n [ nu ] . The stationary solution of (22) can be found in the form of a potential function. Note that it indeed has thesame form as the distribution that we observed earlier (3) to maximize entropy:(23) u ( n ) = 1 Z n exp (cid:26) rn − λn m log( n )) (cid:27) = 1 Z v ( n )e α · A , where v ( n ) = n is the baseline distribution (the stationary solution without any forces acting on the system), A = ( n, − n , log( n )) is a set of observables, and α = ( r, λ, m ) is a set of the ecological forces driving thesystem. The potential function α · A consists of the effects of growth, density regulation, and migration.We assume that migration is strong m > / v ( n ) is not integrable on Ω X =(0 , ∞ ), the function u ( n ) is integrable. The expectations of the observables have biologically meaningfulinterpretations, and can, in principle, be measured. In our case, (cid:104) n (cid:105) corresponds to the expected populationsize, (cid:104) n (cid:105) to the second moment of population size, and the third term, (cid:104) log( n ) (cid:105) is the logarithm of thegeometric mean of the population size. The normalizing constant Z , which is the function of the effectiveforces α , plays an important role, as a generating function for quantities of interest [14](24) ∂ log( Z ) ∂ (2 α j ) = (cid:104) A j ( n ) (cid:105) , ∂ log( Z ) ∂ (2 α i ) = Cov( A i ( n ) , A j ( n )) = C i,j . Given a set of forces α , the system evolves to a stationary distribution (23) that maximizes entropy withconstraints on the observables, where 2 α serve as the Lagrange multipliers. We are interested in how thedynamics change when the set of forces changes in time, and in the most extreme case when the set of initial orces α change rapidly to a new set of values α . The observables will evolve towards the new stationarystate, which creates a path between α and α in the space of effective forces.Under the diffusion approximation we can derive ordinary differential equations (similarly to (5) for thechanges in the mean of the observables A = ( n, − n , log( n ))dd t (cid:104) n (cid:105) = r (cid:104) n (cid:105) − λ (cid:104) n (cid:105) + m, (25) dd t (cid:104) n (cid:105) = (2 m + 1) (cid:104) n (cid:105) + 2 r (cid:104) n (cid:105) − λ (cid:104) n (cid:105) , (26) dd t (cid:104) log( n ) (cid:105) = r − λ (cid:104) n (cid:105) + (cid:18) m − (cid:19) (cid:28) n (cid:29) , (27)where the choice of A follows from the stationary form (23). This dynamical system is not closed, yet, wemay apply the DME method to derive a 3-dimensional approximation for the dynamics of effective forces α ∗ of the form (12) with a particular form of matrices B α ∗ and C α ∗ derived in Appendix B. Note that themethod is fully general and may be applied for arbitrary functions α ( t ), capturing non-stationary ecologicalsituations.4.2. Numerical example.
To understand the relationship between the dynamics of the original system(21) and the dynamics of the reduced system we simulated individual population size trajectories using theEuler-Maruyama method and then compared them to the predictions of the DME method. In figure 3A,we see the simulated trajectories for three sets of initial conditions. The initial conditions were not fixed,instead, they were randomly drawn from a stationary initial distribution, parametrized by the growth rate( r ), the strength of density regulation ( λ ), and migration ( m ). Starting in equilibrium, we changed theenvironmental forces abruptly to new values α = ( r, λ, m ) at time t = 0. This forced the system out ofequilibrium and shifted the trajectories toward the new equilibrium, independent of the initial condition.Instead of using the original stochastic differential equation, in the DME we follow the dynamics of theeffective forces α ∗ = ( r ∗ , λ ∗ , m ∗ ) as they move along to the new equilibrium shown in figure 3B. Note, thatthe parameter space is 3 dimensional, but only a 2 dimensional projection is presented here. A single pointin this space describes the full distribution of population sizes.Figure 3C shows that dynamics of the effective forces in the DME approximation of the stochastic logisticmodel with migration is irreversible, i.e., the trajectory in the space of effective forces when the systemchanges from α to α is different from changing α to α . In both cases the system was initialized withthe stationary distribution and the forces were changed at time t = 0.But how close is the distribution approximated by DME to the real distribution? Despite the simple formof our equations, it is not possible to solve it explicitly analytically. Thus, we compared the numericallycomputed distributions for the original (i.e., exact) problem with the numerically computed distributionsobtained by the DME approximation (in figure 3D). We used three approximations to solve the originalproblem: (1) using the trajectories from figure 3A (approaching the exact distribution when the time stepis small and the number of trajectories is large), (2) numerically solving the FPE for the process by thenative solver of Mathematica , and (3) using the transition matrix method, where we follow a Markov-chain,a continuous time birth-death process on the discrete space of non-negative integers. On the other hand, weused an Euler scheme to solve the DME system.Although all methods are approximate, the original problem can be solved with any precision usingmethods (1-2) and thus we may focus here on the accuracy of the DME method itself. We find that theyare in a good agreement with each other, the only exception being the transition matrix method, which isdefined on a discrete space (and to retain biological meaning also has a slightly different variance from (21)).We compared the distributions at different time points in panel D. We observe that although the transitionin the observable quantities is rather slow and monotonic, the changes in the effective forces can be abruptand non-monotonic. Note, that the DME method does not guarantee that the microscopic population sizedistributions are identical, in fact, they can differ substantially. Nevertheless, the DME method aims tocapture the agreements between the macroscopic variables. Figure 3 (E) shows all three key observables andshows that the DME method is in an excellent agreement with the model.
10 20 30 40020406080100120 t n ( A ) Stochastic simulation - - - - r λ ( B ) DME dynamics of effective forces m ( C ) Irreversibility Ψ ( n ) ( D ) Comparison of distributions t = = = = ( discrete ) t m o m en t s o f n ( E ) Time evolution of observables n n / [ n ] MaxEnt solutionFPE solution individual trajectoriesTM solution ( discrete ) Figure 3. (A) Numerical simulations of stochastic population dynamics on a sin-gle island with immigration. Parameters are α = { r, λ, m } = { . , . , } . Weused initial conditions, α , = { . , . , } (black), α , = { . , . , } (blue), and α , = { . , . , } (green). (B) Corresponding dynamics of the effective forces projectedto the ( r, λ ) space. (C) Irreversibility of the process: 2D projections of the trajectories be-tween α = { . , . , } and α , = { . , . , } and reversed are not the same. (D)Histograms of population sizes at t = 1 , , ,
40 with initial condition α , = { . , . , } (black curves in panels A-B). The numerical solution of the corresponding FPE, the dis-crete transition matrix prediction, and the DME all show a close match. (E) The threeobservables n , log( n ), and n /
2. Code in the Electronic Supplementary Information.In the previous example we demonstrated that the method works well in the most extreme case, when theforces in the dynamics change abruptly. This is surprising, as the DME approximation is based on a quasi-stationary assumption, which would intuitively work only when the forces are changing slowly. This maysuggest that the method will perform even better under slow environmental changes, which in reality maybe more likely to happen than an abrupt change. Temporal differences in the environment can be abrupt,causing populations to become maladapted and possibly drive them to extinction. Less rapid environmentalchanges can be observed on various timescales, for example the warming of the oceans [28] or the yearly cycleof seasons [29], resulting in different migration rates throughout the year due to varying resource abundance.In figure 4 we compare periodic changes for three scenarios: an abrupt, a slow and a fast but continuouschange. We found again that the solution of DME is in a good agreement with that of the FPE, with thetemporal dynamics in the DME and the moments in the FPE equation being indistinguishable by eye. Wealso see that the solution of the non-equilibrium dynamics lags behind the equilibrium of the environment,and that the amount of this lag depends on the speed of change of the ecological forces. Moreover, fasterenvironmental oscillations also result in a smaller range of effective forces. In the extreme case of very fastenvironmental changes (e.g. oscillations with large frequency, or environmental temporal noise of a fixed PE solutionMaxEnt solutionoptimum t i m e t i m e t i m e (A) Abrupt changes (B) Smooth slow changes (C) Smooth fast changes Figure 4.
Periodic changes in the carrying capacity between 20 and 50. The system startsfrom equilibrium with parameters { . , . , } (as in Figure 3), then periodic shifts occurbetween { . , . , } and { . , . , } . The equilibrium distribution of population sizeis shown as it changes in time (background colors). The black dashed line is the expectationof these distributions, the black solid line shows the solution of the DME, whereas the whiteis the solution of the FPE.variation) one expects that effective forces will stay almost stationary as the convergence to equilibrium ismuch slower than the timescale of environmental fluctuations.5. Discussion
We presented an application of a DME method, which helps reduce complexity of the stochastic processby linking the microscopic quantities to the macroscopic observables using a dynamic modification of themaximum entropy method.We first studied the OU process after a rapid change in its parameters. This example demonstrates ourmain strategy for understanding non-stationary dynamics – instead of following the full stochastic dynamicsto follow just the key observables (in this case the first two moments), which change deterministically. Theobservables and the forces acting on them were identified from the potential form of the stationary solution,which is also a solution of a maximum entropy problem and for the OU process has a Gaussian form. Thedynamic problem was solved by the DME method that uses the stationary ansatz but allows the forces tochange in time to best approximate the dynamics of observables. We derived a two-dimensional dynamicalsystem of the effective forces that characterize the solution of the OU dynamics. Despite the intricacies ofthe DME approximation, the DME dynamics coincides with an exact solution of the OU process. This isbecause the dynamical equations for the first two moments are closed and thus the solution is Gaussian atall times. However, even though the OU dynamics are linear, the effective forces solve a nonlinear system ofordinary differential equations.The focus of our work is on the stochastic island model, represented by nonlinear population dynamicsbased on a logistic equation supplemented by migration from other islands. The key parameters of theproblem, the intrinsic growth rate, the carrying capacity, and the migration rate, are in general functions oftime, reflecting temporal environmental changes, which make the problem out-of-equilibrium. Nonlinearity ofthe process results in the dynamics of moments (cid:104) n k (cid:105) which are not closed for any k. Therefore we used DMEto derive a 3-dimensional nonlinear dynamical system for the effective forces using DME. The associatedobservables are no longer just the first three moments but include (cid:104) n (cid:105) , (cid:104)− n / (cid:105) and (cid:104) log n (cid:105) . Unlike for theOU model, the DME approximation of the island model with migration is no longer exact. The system is notfully explicit but contains terms using hypergeometric functions. Nevertheless, it can be solved for dynamicenvironmental forces using standard numerical solvers.We found that the effective forces in the DME approximation lag behind the true environmental forces,which is more pronounced when the environmental forces change faster. However, even in cases of rapid hanges of the environmental forces the observables in the DME approximation are still extremely accurate atall times and thus the effective forces serve as a proxy for the dynamics. When the environmental forces settleto constant values, the effective forces also converge to these values. The DME serves as a change of opticswhere we represent the non-equilibrium dynamics using a series of equilibria parametrized by dynamicaleffective forces. The strength of the DME method in our view is not in speeding up the numerical methodfor solving the problem but in understanding the underlying dynamics in an appropriate low-dimensionalspace.Although it is possible to use the quasi-stationary approximation without the connection to the ME,physics provides us with useful information. By formulating the suitable ME problem for the stationarydistribution we learned which macroscopic quantities are important for the low-dimensional projection ofthe dynamics. Moreover, the observables and the environmental forces are dual variables in the appropriateME formulation. While the observables enter the variational problem via the constraints, the forces are thecorresponding Lagrange multipliers in ME.In addition, the DME method can be placed into the arsenal of methods for non-equilibrium dynamics, aswe have shown in figure 1. It is built from the stationary ME method but unlike the MC method, which isessentially a ME method applied on temporal trajectories, it uses the FPE to establish relationships betweenthe time points.One of the most striking properties of the DME method is its accuracy on the macroscopic level. This issurprising because the quasi-stationary assumption suggests validity of the approach when the applied forceschange adiabatically. However, even for the fast changing forces the approximation stays very accurate. Theunusual form of the DME approximation makes the analytical study of the accuracy of the method a difficultmathematical problem which remains an open problem to this date despite insight provided in [21].This work outlines the first step towards studying more complex questions where the complexity ofthe problem is prohibitive for studying the full problem. In particular, our future goal is to explore eco-evolutionary dynamics where the ecological and population genetic timescales interact. Such interactionhas been studied in [26] but only in the stationary case. Although the existence of an explicit stationarydistribution in principle allows us to explore the dynamics in the non-stationary environment using DME,the structure of the problem poses multiple difficulties that need to be resolved first.The approach may also be suited to stochastic problems in different disciplines. The method is based onthe structure of the problem, in which the stochastic dynamics are described by the FPE and the stationarysolution is explicit. This includes a wide range of problems accross disciplines, for example stochasticcoagulation-fragmentation dynamics when the rates satisfy a detailed balance condition (existence of anexplicit stationary distribution for this problem was shown in [30]). Funding
This work has been supported by the Scientific Grant Agency of the Slovak Republic under the GrantsNos. 1/0755/19 and 1/0521/20.
Appendix A. DynMaxEnt for the OU process, derivation of the B and C matrices We briefly outline the key steps in the derivation of DME. The forces and observables are α = (2 µβ, − β )and A = ( x, x ). The expectations (cid:104) A i (cid:105) follow a closed system of ODEs(A.1) (cid:104) A (cid:105) (cid:48) = (cid:18) / (cid:104) x (cid:105)(cid:104) x (cid:105) (cid:104) x (cid:105) (cid:19) α + (cid:18) σ (cid:19) = B α + V . To apply the DME approximation we assume that at every time there are effective forces α ∗ such that B α ∗ α ∗ + V ∗ = 0 (with the moments in matrix B evaluated at the stationary distribution ¯ u α ∗ ). Usingstationarity condition V ∗ = − B α ∗ α ∗ to substitute V by V ∗ and B by B α ∗ in (A.1) we obtain(A.2) dd t (cid:104) A (cid:105) = B α ∗ ( α − α ∗ ) . The matrix B α ∗ can be expressed in terms of the effective forces α ∗ , although we will write most of theexpressions in terms of µ ∗ and β ∗ (the transformation from α ∗ to ( µ ∗ , β ∗ ) is regular). The matrix B α ∗ can e expressed as(A.3) B α ∗ = (cid:18) / µ ∗ µ ∗ µ ∗ ) + σ /β ∗ (cid:19) . Finally, we change variables in (A.2) using the scaled covariance matrix { C α ∗ } ij = Cov( A i , A j ) /σ to obtaindynamics of α ∗ (A.4) C α ∗ = d (cid:104) A (cid:105) α ∗ d α ∗ = 12 β (cid:18) µ µ σ /β + 4 µ (cid:19) , C − α ∗ = 2 βσ (cid:18) βµ + σ − µβ − µβ β (cid:19) , where for simplicity we dropped the ∗ notation in the above equations (all coefficients µ and β are understoodas µ ∗ and β ∗ ). Plugging this into (A.2) (this time keeping all ∗ symbols) leads tod α ∗ d t = C − α ∗ B α ∗ ( α − α ∗ ) = 2 β ∗ (cid:18) − µ ∗ (cid:19) (cid:18) βµ − β ∗ µ ∗ β ∗ − β (cid:19) . (A.5)After some algebraic manipulation the dynamics of α ∗ becomesdd t (2 µ ∗ β ∗ ) = 2 β ∗ ( µβ − µ ∗ β ∗ ) + 2 µ ∗ β ∗ ( β − β ∗ ) , (A.6) dd t ( − β ∗ ) = 2 β ∗ ( β ∗ − β ) , (A.7)Finally we transform the dynamics of α ∗ to µ ∗ , β ∗ d µ ∗ d t = β ( µ − µ ∗ ) , (A.8) d β ∗ d t = 2 β ∗ ( β − β ∗ ) . (A.9)This system is identical to (19)-(20). Appendix B. DynMaxEnt for the island model, derivation of the B and C matrices Equation 25 can be written using the matrix notation: (cid:104) A (cid:105) (cid:48) = (cid:104) n (cid:105) (cid:104) n (cid:105) (cid:104) n (cid:105) (cid:104) n (cid:105) (cid:104) n (cid:105) (cid:104) n (cid:105) (cid:10) n (cid:11) α + (cid:104) n (cid:105)− (cid:10) n (cid:11) = B α + V . (B.1)At each time point, we approximate the elements of B and V using the stationary approximation B α ∗ α ∗ + V ∗ = 0. Substituting V = − B α ∗ α ∗ into (B.1) and using that B ≈ B α ∗ , we obtain ∂ (cid:104) A i ( n ) (cid:105) ∂t ≈ (cid:80) j B ∗ i,j ( α j − α ∗ j ). Change of variables yields(B.2) d α ∗ d t = (cid:20) d (cid:104) A (cid:105) α ∗ d α ∗ (cid:21) − d (cid:104) A (cid:105) α ∗ d t . The expectations of various functions of variable n appearing in matrices B α ∗ and C α ∗ can be expressedanalytically under the condition that the migration rate is not too low ( m > ). Let us call the k th momentof the stationary distribution G ( k ), this can be expressed analytically in terms of hypergeometric functions(suppressing ∗ notation) G ( k ) = (cid:90) ∞ n k v ( n )e α A = (cid:90) ∞ n k +2 m − exp {− λn + 2 rn } = 12 λ ( − − k − m ) × (B.3) × (cid:18) √ λ Γ (cid:18) k m (cid:19) F (cid:18) k m, , r λ (cid:19) + 2 r Γ (cid:18) k + 12 + m (cid:19) F (cid:18) k + 12 + m ; 32 ; r λ (cid:19)(cid:19) , (B.4)if Re( k + 2 m/γ ) >
0. Using the function G , all the moments of interest can be expressed as(B.5) G (0) = Z , G (1) G (0) = (cid:104) n (cid:105) , G (2) G (0) = (cid:104) n (cid:105) , G (3) G (0) = (cid:104) n (cid:105) , G ( − G (0) = (cid:28) n (cid:29) . hus B α ∗ = 1 G (0) G (1) − G (2) 1 − G (2) G (3) G (1)1 − G (1) G ( − . (B.6)Furthermore, we can express (cid:104) log( n ) (cid:105) analytically by taking the j th derivative of G ( k ) with respect to m :(B.7) H ( k, j ) = E ( n k log( n ) j ) = 12 j ∂ ( j ) G k ∂m j . The covariance matrix of the observables evaluated at the quasi-stationary distribution parametrized by theeffective forces can then be written as(B.8) C α ∗ = G (2) G (0) − G (1) G (0) (cid:16) G (1) G (2) G (0) − G (3) G (0) (cid:17) H (1 , H (0 , − G (1) H (0 , G (0) H (0 , (cid:16) G (1) G (2) G (0) − G (3) G (0) (cid:17) (cid:16) G (4) G (0) − G (2) G (0) (cid:17) (cid:16) G (2) H (0 , G (0) H (0 , − H (2 , H (0 , (cid:17) H (1 , H (0 , − G (1) H (0 , G (0) H (0 ,
0) 12 (cid:16) G (2) H (0 , G (0) H (0 , − H (2 , H (0 , (cid:17) H (0 , H (0 , − H (0 , H (0 , . References [1] Jaynes ET. 1957 Information theory and statistical mechanics. Phys. Rev. : 620.[2] Bialek W, Cavagna A, Giardina I, Mora T, Silvestri E, Viale M, Walczak AM. 2012 Statistical me-chanics for natural flocks of birds. Proc. Natl. Acad. Sci. U.S.A. : 4786–4791.[3] Schneidman EM, Berry J, Segev R, Bialek W. 2006 Weak pairwise correlations imply strongly correlatednetwork states in a neural population. Nature : 1007–1012.[4] Mora T, Walczak AM, Bialek W, Callan CG. 2010 Maximum entropy models for antibody diversity.Proc. Natl. Acad. Sci. U.S.A. : 5405–5410.[5] Thompson JN. 1998 Rapid evolution as an ecological process. Trends in Ecology & Evolution, (8),329-–332.[6] Schoener TW. 2011 The newest synthesis: understanding the interplay of evolutionary and ecologicaldynamics. Science, (6016), 426–429.[7] Uecker H, Otto SP, Hermisson J. 2014 Evolutionary rescue in structured populations. The AmericanNaturalist, (1), E17–E35.[8] Sachdeva H. 2019 Effect of partial selfing and polygenic selection on establishment in a new habitat.Evolution, (9), 172–1745.[9] Polechova J. 2018 Is the sky the limit? on the expansion threshold of a species’ range. PLoS Biology, (6), e2005372.[10] Jaynes ET. 1980 The minimum entropy production principle, Annu. Rev. Phys. Chem. , 579.[11] Press´e S, Ghosh K, Lee J, Dill KA. 2013 Principles of maximum entropy and maximum caliber instatistical physics. Rev. Mod. Phys. : 1115.[12] Dixit PD et al. . 2018 Perspective: Maximum caliber is a general variational principle for dynamicalsystems. J Chem. Phys. : 010901.[13] Ghosh K et al. :213–38.[14] Barton NH, de Vladar HP. 2009 Statistical mechanics and the evolution of polygenic quantitativetraits. Genetics : 997–1011.[15] de Vladar HP, Barton NH. 2011 The statistical mechanics of a polygenic character under stabilizingselection, mutation and drift. J. R. Soc. Interface (58): 720–739.[16] Bod’ov´a K, Tkaˇcik G, Barton NH. 2016 A General Approximation for the Dynamics of QuantitativeTraits. Genetics (4): 1523–1548.[17] Hick P, Stevens G. 1987 Approximate solutions to the cosmic ray transport equation the maximumentropy method. Astron. Astrophys. : 350–358.[18] Tsallis C, Bukman DJ. 1996 Anomalous diffusion in the presence of external forces: Exact time-dependent solutions and their thermostatistical basis Phys. Rev. E : R2197.[19] Plastino AR, Miller H, Plastino A. 1997 Minimum kullback entropy approach to the fokker-planckequation. Phys. Rev. E : 3927.
20] Plastino AR, Plastino A. 1997 Statistical treatment of autonomous systems with divergenceless flowsPhysica A : 458–476.[21] Bod’ov´a K, Haskovec J, Markowich P. 2018 Well posedness and maximum entropy approximation forthe dynamics of quantitative traits. Physica D (4): 108–120.[22] Segel LA, Slemrod M. 1989 The quasi-steady-state assumption: a case study in perturbation. SIAMRev. : 446–477.[23] Goeke A, Walcher S. 2013 Quasi-steady state: Searching for and utilizing small parameters. In RecentTrends in Dynamical Systems , pp. 153–178, Springer.[24] Koll´ar R, ˇSiˇskov´a K. 2015 Extension and justification of quasi-steady-state approximation for reversiblebimolecular binding. Bull. Math. Biol. : 1401–1436.[25] De Vladar HP, Barton NH. 2011 The contribution of statistical physics to evolutionary biology. Trendsin ecology & evolution (8), 424-432.[26] Sz´ep E, Sachdeva H, Barton NH. 2020 Polygenic local adaptation in metapopulations: a stochasticeco-evolutionary model. bioRxiv. doi: https://doi.org/10.1101/2020.06.16.154245[27] Ewens WJ. 2012 Mathematical population genetics 1: theoretical introduction. Springer Science &Business Media, Vol. 27, Chap. 4-5.[28] Munday PL, Warner RR, Monro K, Pandolfi JM, Marshall DJ. 2013 Predicting evolutionary responsesto climate change in the sea. Ecology letters, (12), 1488–1500.[29] Kingsolver JG, Buckley LB. 2017 Evolution of plasticity and adaptive responses to climate changealong climate gradients. Proceedings of the Royal Society B: Biological Sciences, (1860), 20170386.[30] Durret R, Granovsky BL, Gueron S. 1999 The equilibrium behavior of reversible coagulation-fragmentation processes. J. Theor. Prob. (2): 447–474.(2): 447–474.