Dynamic maximum entropy provides accurate approximation of structured population dynamics

Realistic models of biological processes typically involve interacting components on multiple scales, driven by changing environment and inherent stochasticity. Such models are often analytically and numerically intractable. We revisit a dynamic maximum entropy method that combines a static maximum entropy with a quasi-stationary approximation. This allows us to reduce stochastic non-equilibrium dynamics expressed by the Fokker-Planck equation to a simpler low-dimensional deterministic dynamics, without the need to track microscopic details. Although the method has been previously applied to a few (rather complicated) applications in population genetics, our main goal here is to explain and to better understand how the method works. We demonstrate the usefulness of the method for two widely studied stochastic problems, highlighting its accuracy in capturing important macroscopic quantities even in rapidly changing non-stationary conditions. For the Ornstein-Uhlenbeck process, the method recovers the exact dynamics whilst for a stochastic island model with migration from other habitats, the approximation retains high macroscopic accuracy under a wide range of scenarios in a dynamic environment.


Introduction
Conceptual understanding of realistic problems in applied sciences is often hindered by the curse of complexity, with quantities of interest coupling to finer features. Due to their multiscale character even simple questions lead to exploration of the full complexity of the system. But how can we ever understand the processes around us if incremental learning is impossible?
Statistical mechanics provides a clever way to understand complex multiscale problems by linking processes on 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 Fig 1A. The method gained 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. Adaptation of populations to spatial and temporal ecological gradients is an example of a complex non-equilibrium processes in which ecological and evolutionary processes interact [5][6][7][8][9]. How can a simple concept of ME be applied to such systems? The most straightforward approach is to use ME on dynamic trajectories, forcing constraints on the dynamical features ( Fig 1B). This approach, called maximum caliber (MC), introduced by [10] is suitable for 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 microscopic distribution, given knowledge of a few key macroscopic observables. (B) MC is analogous to ME, however, each data point represents a trajectory x(t). MC connects the microscopic distribution 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, given by the FPE, which reduces the full problem to a low-dimensional dynamics. This reduction is a consequence of a ME ansatz; the approximation at each time � u DME ðαðtÞÞ solves the ME problem (stationary form in the FPE), where the dynamics of the effective forces α are systematically derived from the FPE. inference of models and their parameters from dynamic data. Comprehensive reviews of MC [11][12][13] provide multiple examples where the method has been successfully applied to non-stationary biological processes.
While such an inverse approach is useful for understanding temporal data, in many cases data are not available but the dynamics, although often extremely complicated, are known to be accurately described by the Fokker-Planck equation (FPE). The aim of our work is to demonstrate usefulness of a theoretical dimensionality-reduction technique, which approximates possibly many-dimensional stochastic dynamics to a low-dimensional deterministic dynamics of a few key observable quantities (Fig 1C). The approximation combines a static ME ansatz in the FPE equation with a quasi-stationary assumption. We refer to this method as dynamic ME (DME). DME has been used before to solve problems in quantitative genetics [14][15][16] and independently in cosmology [17][18][19][20]. While the problems in quantitative genetics are formulated using a linear FPE and ME uses relative entropy, the applications in physics are formulated using nonlinear FPE and the ME is based on Tsallis entropy. This work studies linear dynamics, characterized by a linear FPE.
The most surprising feature of DME is its accuracy. The method, although derived from the assumption of quasi-stationarity, remains extremely accurate even far-from-equilibrium. Although numerical evidence of the method's accuracy has been demonstrated in [14][15][16]21], explicit analysis is a challenging mathematical problem. This is analogous to problems in reaction kinetics, where reduction of dynamics based on the quasi-stationary assumption is often used but its validity has been shown rigorously only for a few basic systems using singular perturbation [22][23][24]. In comparison, the accuracy of DME still remains a mystery.
We will address this puzzle by focusing on two processes: the Ornstein-Uhlenbeck (OU) process, and a logistic model of population growth in a continent-island model. The OU process is studied as an example of a dynamics, for which DME gives an exact solution even in a fully general (non-stationary) form, where the two key parameters are functions of time. We derive the exact solution as well as the DME approximation and show that the dynamics are the same. The stochastic island model, as one of the simplest, yet nonlinear, models in population dynamics is studied as a first step to understand more complex population models where ecological processes interact with evolutionary processes. Despite its simple form the DME is no longer exact, and even fails when migration to the island is less than a threshold 1/2. On the other hand, using numerical simulations we show that for migration rate above this threshold the approximation is extremely accurate even far from equilibrium.

Dynamic maximum entropy
Here we present a dynamic maximum entropy (DME) method to approximate stochastic dynamics by a FPE [14][15][16][17][18][19][20][21]. The method is based on a combination of ME in statistical physics [1], which solves the stationary 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 stationary distribution, 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). DME was introduced in population genetics to understand how quantitative traits change in time in the presence of various evolutionary mechanisms without resolving details about the dynamics of the underlying gene frequencies. However, the origins reach back to [25] studying genetic algorithms for finding the low-energy state in the Ising spin glass by an analogy with statistical mechanics. The problem reduced to an infinite-dimensional dynamics of cumulants, which were solved numerically for the truncated system. Later, the authors of [26] derived the cumulant dynamics for the population genetics problem with population under selection, mutation and drift. To close the infinite-dimensional dynamical system the authors used perturbation theory in the weak selection regime and separation of time scales in the weak mutation regime. The problem of obtaining closed macroscopic dynamics was systematically resolved in [14] by introducing the DME approximation for the dynamics of quantitative traits under directional selection, mutation and drift. The low-dimensional macroscopic dynamics were obtained by considering suitable macroscopic variables, consistent with the ME solution of the stationary problem, and by using the local equilibrium approximation. The method was further applied to a polygenic trait under stabilizing selection, mutation and drift in [27] and its limitations in the small-mutation regime were resolved in [16] by distinguishing the bulk and the boundary layers of the microscopic distributions of allele frequencies.
Although accuracy of the method is still not rigorously understood, the authors of [21] showed an exponential convergence to equilibrium in the FPE by finding a positive spectral gap and investigated an alternative formulation of the DME, which is not exact in the stationary case, but leads to a lower error in the dynamical situation and avoids the problems with small mutation rate.
Independent use of the method in statistical physics focused on exact and approximate solutions of a nonlinear FPE arising in cosmology and in other areas of physics. We traced the first relevant connection to [28], where a relationship was established between the linear Fokker-Planck equation and maximum entropy problem, which uses the generalized Tsallis entropy by finding a family of linear FPEs, whose stationary solutions solve the ME problem with the Tsallis entropy. The use of the generalized entropy led to further studies of problems, mathematically described by nonlinear Fokker-Planck equations. Namely, the authors of [20] found particular time-dependent solutions of the nonlinear FPE with linear drift (generalization of the Ornstein-Uhlenbeck process) and power-law noise, which are solutions of the ME problem with the Tsallis entropy. Their ideas were further developed in [29], introducing an approximate ME approach for the study of the nonlinear FPE, based on the Tsallis entropy, called nonextensive maximum entropy approach. This approximation is analogous to the DME approximation, although the class of problems, as well as the definition of the entropy differs.
Although in most of the previous work the derivation of the DME method was presented phenomenologically, a rigorous derivation of the approximation can be found in [21]. Here we provide a comprehensive summary of the method.
Assume stochastic dynamics with a potential in the Langevin form , and x(0) = x 0 . The potential in the first term (in the brackets) is a linear combination of time-dependent forces α i (t), 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][15][16]27]) focused on examples in population genetics where x k corresponds to the frequency of a certain gene, affecting some quantitative trait. This frequency depends on evolutionary processes, e.g. selection, mutation, and inherent stochastic fluctuations, described by the forces α i (t). The distribution u(t, x) follows dynamics described by the FPE @uðt; xÞ @t ¼ À which can be also expressed in the flux form @ t uðt; xÞ ¼ À P N k¼1 @ x k J k ½t; x�. This FPE is complemented with no-flux boundary conditions J k [t, x] = 0 at x k 2 @O X and the initial condition u(0, x) = u 0 (x).

Stationary solution and ME
If the vector of forces αðtÞ ¼ ða 1 ðtÞ; . . . ; a d ðtÞÞ 2 R d is time-independent then at large times the dynamics approach a stationary distribution, parametrized by the vector of limiting forces with the normalization coefficient (i.e., the partition function) where A = (A 1 , . . ., A d ) is a vector function of the state variables x which the forces α i act on.
Using the terminology of statistical physics we refer to functions A i as observables, as their expectations in problems in 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 problem with 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 α. The dynamics of the expectation hA j i follows from the FPE. Using notation This forms a system of ordinary differential equations for hAi, which is generally not closed due to nonlinearity of the functions A j (x). Next we define a logarithmic relative entropy where u(t, x) is a solution of the FPE at time t � 0 and x k 2 O X . For any t � 0 relative entropy in Eq 6 has a maximum at α = α � , which can be obtained by solving a set of first-order conditions The above intriguing relationship states that if for a given time t there exists a maximum of the relative entropy in Eq 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 time. This simple fact, shown previously in [21] and in a slightly different form in [14][15][16] suggests that instead of the full representation of the problem using FPE one could trace only the d-dimensional dynamics of α � , that parametrize the approximate solution by the form in Eq 3. Furthermore, the two representations agree in terms of the expectations hA i i at a given time. Concavity of the relative entropy is implied by the following relationship analogous to similar expressions in [14][15][16]21] stating that the Hessian of the relative entropy is positive semidefinite. Note that the solvability of the Eq 7 (for α � ) follows from Eq 8 when the covariance matrix is positive definite [21] (it is always positive semidefinite). Note that the DME can be also defined using the relative entropy, which compares the distribution u with the reference distribution ϕ, where ϕ is a solution of the problem in the absence of forces α. This approach was used in the population genetic applications where ϕ represented the neutral distribution in the absence of selection and mutation [14,16,27]. ME is then formulated as a constrained optimization where each force α i (evolutionary/ecological) enters the problem in the form of the Lagrange multiplier, corresponding to a constraint on a specific complementary macroscopic quantity A i . The alternative formulation of the ME leads to the same mathematical outcome as the one presented in this section. Moreover, it allows us to keep some of the forces constant throughout the evolutionary timescales considered by including them in the reference distribution, while others, which are dynamically adapting, added through the constraints.

Dynamical approximation
We have established a relationship between the solution of the full stochastic dynamics in Eq 2 and a stationary form u α � parametrized by suitable effective forces α � following ME. However, as we demonstrated in Fig 1A, ME is applicable only to static problems. When the system is out-of-equilibrium, we need to establish a dynamic relationship between the values α � (t 1 ) and α � (t 2 ) for t 1 6 ¼ t 2 by using the information captured by the FPE.
To derive the DME approximation of Eq 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 in Eq 5 become where h�i u represents expectation through distribution u and hf ðxÞi α ≔ hf ðxÞi � u αðtÞ . Now we make two key assumptions. First, we assume that the residual terms in the bracket of Eq 9 are small and we neglect them. In addition, we also impose a quasi-stationarity approximation, assuming that α � are chosen to satisfy the equilibrium relationship for all t > 0. Both steps are easy to justify if the forces α(t) are slowly changing (in the adiabatic regime) and the solution of the FPE is thus close to an equilibrium form in Eq 3. However, its validity when out-of-equilibrium is not clear. We use Eq 10 to replace V α by V α � and B α by B α � in Eq 9 to approximate Finally, to obtain a closed dynamical system for α � we use the chain rule, noting that

@t . Combination of Eqs 7 and 8 imply that differentiation of an expectation
together with a parametric form in Eq 3 represent DME approximation of dynamics in Eq 2. Solution of Eq 12 can be plugged into the stationary parametric form in Eq 3, which allows us to not only study the accuracy of the key moments used in DME, but also to compute any statistical feature of the approximate solution and compare it with the exact solution.
Note that Eq 12 can be solved for any prescribed continuous function α(t) and in the most extreme case, the forces α(t) can contain step changes, which clearly violate the quasi-stationarity assumption. However, all previously studied applications of DME approach showed that the approximation captures extremely well the expectations of the key functions even when the forces change rapidly. This is one of the most 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) typically involve high-dimensional stochastic dynamics with nonlinearities and coupling terms [14][15][16]. However, for simplicity we consider examples leading to a simpler one-dimensional form with a remark that more complex models can be analyzed using this approach as long as the stationary distribution of the FPE is explicit. However, even when x is a scalar A and α are vectors in all problems studied here: these vectors summarise the infinite-dimensional distribution of x.

The model
The importance of examples where the DME becomes exact may lead us to the understanding of why the approximation works even for more general cases. In the area of nonlinear FPE such an example was provided by [20] who analyzed a stochastic process with a linear advection and power-law noise, i,e, a generalization of the Ornstein-Uhlenbeck process. They showed that the nonextensive maximum entropy method is exact for the studied process.
Here we outline a simple example of stochastic dynamics where DME reproduces the exact dynamics, namely the standard OU process with linear relaxation to an equilibrium in the presence of constant Gaussian noise (examples include a particle under friction, animal motion, financial time series, etc.). The presentation in this section is partly pedagogical, although we are not aware of other works where the OU process is studied within the context of DME. The state variable x in the OU process dynamically adapts to the value μ at a speed given by β, further affected by a white noise of a magnitude σ We consider β(t) and μ(t) to be time-dependent, with long-time limits β 1 and μ 1 and unless otherwise stated we consider σ(t) = σ 0 to be constant. Note that σ(t) can be replaced by a constant using transformation of time in Eq 13 (when σ(t) is smooth and invertible) without changing the character of the problem. The stationary distribution of the Eq 13 can be obtained from the FPE, which describes time-evolution of the probability distribution of x, denoted by u(t, x) by letting t ! 1. Defining this yields where Z ¼ expðm 2 1 =s 2 0 Þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 2ps 2 0 =b 1 p is the normalization factor, α 1 = (2μ 1 β 1 , −β 1 ) and A = (x, x 2 ). If the initial condition in Eq 13 is the stationary distribution corresponding to the parameters (β 0 , μ 0 , σ 0 ), i.e., as a Gaussian x 0 � N ðm 0 ; s 2 0 =2b 0 Þ, then the solution of Eq 14 is a Gaussian N ðmðtÞ; vðtÞÞ for every t � 0. This well-known result can be found for instance in [30,31]. The general form of this solution is u( where the functions f, g, h satisfy the following dynamical system Note that the last equation is just the deterministic version of the OU process. In the special case when the coefficients β(t) = β 1 and μ(t) = μ 1 for all t > 0 while β 1 6 ¼ β 0 , μ 1 6 ¼ μ 0 and σ 0 are constants the mean m(t) and the variance v(t) solve the system of ODEs with initial conditions m(0) = μ 0 and vð0Þ ¼ s 2 0 2b 0 . This system has an explicit solution which satisfies m(t)!μ 1 and vðtÞ ! s 2 0 =2b 1 as t ! 1. The stationary solution � u α 1 ðxÞ in Eq 16 solves a variational ME problem with the relative entropy defined in Eq 6 using forces α 1 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 in Eq 16, due to the fact that time-dependence in σ(t) can be removed by transforming time in the dynamics. This allows us to keep the volatility of the process fixed in the DME approach.
We briefly outline the key steps in the derivation of DME. The forces and observables are α(t) = (2μ(t)β(t), −β(t)) and A = (x, x 2 ). In the following calculations we assume temporal dependence of α but we suppress its notation. The expectations hA i i follow a closed system of ODEs To apply the DME approximation we assume that at every time there are effective forces α � = (2μ � β � , −β � ) such that B α � α � + V � = 0 (with the moments in matrix B α � evaluated at the stationary distribution � u α � ). Using the stationarity condition V � = −B α � α � to substitute V by V � and B by B α � in Eq 23 we obtain In essence, DME approximates the general solution of Eq 14 with the Gaussian initial condition by a stationary form � u α � ðtÞ with time-dependent parameters α � (t). The matrix B α � can be expressed in terms of the effective forces α � , although we will write most of the expressions in terms of μ � and β � (the transformation from α � to (μ � , β � ) is regular). The matrix B α � can be expressed as Finally, we change variables in Eq 24 using the scaled covariance matrix fC α � g ij ¼ Plugging this into Eq 24 leads to After some algebraic manipulation we obtain a 2-dimensional dynamical system for the effective forces α � This coupled dynamical system can be transformed into decoupled dynamics of μ � and β � of the form First, note that the dynamical system Eqs 31 and 32 is independent of the noise magnitude σ 0 , which is considered constant. This is expected and follows from the same scaling of the variance of the OU process (exact or DME approximation) and s 2 0 . For the exact OU process, this scaling can be confirmed for instance by the dynamics of the variance in the form dv dt ¼ À 2bðtÞv þ s 2 0 derived from Eq 14. On the other hand, the DME solution has an explicit Gaussian form � u α � ðtÞ , introduced in Eq 15, which leads to the same scaling between its variance and s 2 0 provided β � (t) is independent of σ 0 .
The dynamical system in Eqs 31 and 32 is explicitly solvable since it is decoupled and while the first equation is linear, the second one is logistic. The DME solution is consistent with Eqs 18 and 19 with the relationship between the dynamic variables h(t) = μ � (t) and s 2 0 gðtÞ ¼ b � ðtÞ. The first two moments of the distribution are given by m(t) = μ � (t) and vðtÞ ¼ s 2 0 =ð2b � ðtÞÞ, which are functions of effective forces α � . This is due to the linearity of the OU process, which yields a closed dynamics of the first two moments and thus preserves a Gaussian form of the solution at each time, provided we started with a Gaussian initial condition. Note that the OU process is not the only stochastic process where DME provides an exact solution. As [18] showed, a nonlinear extension of the OU process can be solved exactly using a ME ansatz.  figure legend). Initially, the system is in a stationary state corresponding to parameters μ 0 , β 0 , σ 0 (Gaussian form with parameters x 0 ; s 2 0 =2b 0 ). However, at t = 0 the parameters of the OU process rapidly change, pushing the system out-of-equilibrium. As a response, the distribution of sample trajectories follows a Gaussian form at each time, eventually converging to N ðm 1 ; s 2 0 =2b 1 Þ. In Fig 2B we plot the 2-dimensional DME dynamics of effective forces β � , μ � , which is exact. Each trajectory (the vector field of the dynamical system is plotted as well) represents the complete solution of the FPE for a given parameter choice (at each time it is a Gaussian). The microscopic distribution at four different times in panel C shows an agreement between stochastic simulations (histograms) and microscopic distributions obtained from the DME approach (solid curves). In general, the goal of DME is to approximate the dynamics on the macroscale, thus, we do not expect DME to capture also the microscopic distribution. However, for the OU process DME is exact and thus the method recovers both the macroscale and the microscale properties of the process without loss of precision.

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 other habitats.

The model
In the case of unlimited resources and in the absence of predation, populations would grow indefinitely. However, in natural populations this is not the case: various factors impose bounds on this exponential growth. The logistic growth model describes population size regulation in the absence of demographic stochasticity (i.e., when the population size is infinite). At low population sizes, when resources are abundant and competition is low, the population grows with its intrinsic growth rate, r. However, the total growth rate of the population decreases linearly with increasing population size. In particular, the growth rate is zero when the population is at carrying capacity, K, when each individual replaces itself in each generation. The carrying capacity represents the maximal sustainable population size. We follow the population in a single island, with a migration from other habitats at rate m. In the presence of demographic stochasticity the described population dynamics can be formulated using a stochastic differential equation where n(t) represents the population size at time t, λ = r/K is the density regulation and γ describes the variance in population size. For the sake of simplicity we fix γ = 1 corresponding to a Poisson(1) number of offspring for each individual (with total variance n). Extinction in this stochastic dynamics for m = 0 is unavoidable from a mathematical point of view (as the process is a critical branching process) but can be prevented by migration when m > 0. In general, complex eco-evolutionary interaction requires including changes in population size due to demographic processes, changes occurring in gene frequencies due to selection, as well as the various feedback mechanisms connecting them. Feedback loops between population sizes and gene frequencies can result from migration and hard selection, i.e. when the size of the population depends on its genetic composition. Such questions were studied in [32], but only in a stationary case. Relaxing the assumption of stationarity makes the problem more realistic, and inherently more difficult. The stochastic logistic dynamics with immigration in Eq 33, despite its simplicity, serves as the first step to understanding the biologically more realistic scenario of stochastic models in population dynamics [33]. We are interested in how changes in the environment reflect the dynamics of biological quantities, particularly when the system is out of equilibrium. Since the model is nonlinear, the dynamics of moments, i.e., average population size, etc., are not closed. Nevertheless, the DME method can be applied to reduce the stochastic dynamics to a low-dimensional deterministic dynamics of the key observables.
Based on Eq 33, we find the corresponding FPE describing the time evolution of the probability distribution u(t, n) (which is of the same form as Eq 2): The stationary solution of Eq 34 can be found in the form of a potential function. Note that it indeed has the same form as the distribution that we observed earlier in Eq 3 to maximize entropy: where vðnÞ ¼ 1 n is the baseline distribution (the stationary solution without any forces acting on the system), A ¼ n; À n 2 2 ; logðnÞ À � is a set of observables, and α = (r, λ, m) is a set of the ecological forces driving the system. The potential function α � A consists of the effects of growth, density regulation, and migration. We assume that migration is strong (m > 1/2) so even though the function v(n) is not integrable on O X = (0, 1), the function u(n) is both integrable and bounded for small population sizes n. The expectations of the observables have biologically meaningful interpretations, and can, in principle, be measured. In our case, hni corresponds to the expected population size, hn 2 i to the second moment of population size, and the third term, hlog(n)i is the logarithm of the geometric mean of the population size. The normalizing constant Z, which is the function of the effective forces α, plays an important role, as a generating function for quantities of interest [14] @logðZÞ Given a set of forces α, the system evolves to a stationary distribution in Eq 35 that maximizes entropy with constraints on the observables, where 2α serve as the Lagrange multipliers. We are interested in how the dynamics change when the set of forces changes in time, and in the most extreme case when the set of initial forces α 0 change rapidly to a new set of values α 1 . The observables will evolve towards the new stationary state, which creates a path between α 0 and α 1 in the space of effective forces.
Under the diffusion approximation we can derive ordinary differential equations (similarly to Eq 5 for the changes in the mean of the observables A ¼ n; À n 2 2 ; logðnÞ where the choice of A follows from the stationary form in Eq 35. Eqs 37-39 can be written using the matrix notation: This dynamical system is not closed, yet, we may apply the DME method to derive a 3-dimensional approximation for the dynamics of effective forces α � of the form in Eq 12 with a particular form of matrices B α � and C α � . To do this, we approximate the elements of B and V using the stationary approximation B α � α � + V � = 0. Substituting V = −B α � α � into Eq 40 and using that B � B α � , we obtain @hA i ðnÞi Thus B α � ¼ 1 Gð0Þ Gð1Þ À Gð2Þ Gð0Þ À Gð2Þ Gð3Þ À Gð1Þ Furthermore, we can express hlog(n)i analytically by taking the j th derivative of G(k) with respect to m: The covariance matrix of the observables evaluated at the quasi-stationary distribution parametrized by the effective forces can then be written as The 3-dimensional nonlinear dynamical system in Eq 41 along with Eqs 46 and 48 defines the DME approximation of the island model. Note that the method is fully general and may be applied for arbitrary functions α(t), capturing non-stationary ecological situations.

Failure of the DME method for small migration rate
An apparent failure of the DME for the island model arises in the parameter regime m < 1/2. The stationary form, which is used as a parametric ansatz in the DME, behaves for low population sizes as n 2m−1 leading to a singularity at n = 0 for m < 1/2. Even though the population size distribution is integrable, the average hn −1 i, appearing in the matrix B � is unbounded. A similar problem arises when a quantitative trait is studied under selection, mutation and inherent stochasticity when mutation is weak. In [16] we showed how to resolve this problem by splitting the domain of the independent variable to a bulk and the boundary, where the singularity occures. A similar approach can be applied here as well.

Numerical example
To understand the relationship between the dynamics of the original system in Eq 33 and the dynamics of the reduced system we simulated individual population size trajectories using the Euler-Maruyama method and then compared them to the predictions of the DME method. In Fig 3A, we see the simulated trajectories for three sets of initial conditions. The initial conditions were not fixed, but instead were randomly drawn from a stationary initial distribution, parametrized by the growth rate (r 0 ), the strength of density regulation (λ 0 ), and migration (m 0 ). Starting in equilibrium, we changed the environmental forces abruptly to new values α = (r, λ, m) at time t = 0. This forced the system out of equilibrium 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 the effective forces α � = (r � , λ � , m � ) as they move to the new equilibrium shown in Fig 3B. Note, that the parameter space is 3 dimensional, but only a 2 dimensional projection is presented here. A single point in this space describes the full distribution of population sizes. Fig 3C shows that the dynamics of the effective forces in the DME approximation of the stochastic logistic model with migration are irreversible, i.e., the trajectory in the space of effective forces when the system changes from α 0 to α 1 is different from changing α 1 to α 0 . In both cases the system was initialized with the stationary distribution and the forces were changed at time t = 0.
How close is the distribution approximated by DME to the real distribution? Despite the simple form of our equations, it is not possible to solve it explicitly analytically. Thus, we compared the numerically computed distributions for the original (i.e., exact) problem with the numerically computed distributions obtained by the DME approximation (in Fig 3D). We used three approximations to solve the original problem: (1) using the full stochastic  in panels A-B). The numerical solution of the corresponding FPE, the discrete transition matrix prediction, and the DME all show a close match. (E) The three observables n, log(n), and n 2 /2. Code in S1 Code.
https://doi.org/10.1371/journal.pcbi.1009661.g003 simulation (Fig 3A, distribution at each time approaches the exact distribution when the time step in the stochastic simulation is small and the number of trajectories is large), (2) numerically solving the FPE for the process by the native 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, we used an Euler scheme to solve the DME system.
Although all methods are approximate, the original problem can be solved with any precision using methods (1-2) and thus we may focus here on the accuracy of the DME method itself. We find that they are in a good agreement with each other, the only exception being the transition matrix method, which is defined on a discrete space (and to retain biological meaning also has a slightly different variance from Eq 33). We compared the distributions at different time points in panel D. We observe that although the transition in the observable quantities is rather slow and monotonic, the changes in the effective forces can be abrupt and non-monotonic. Note, that the DME method does not guarantee that the microscopic population size distributions are identical, in fact, they can differ substantially. Nevertheless, the DME method aims to capture the agreements between the macroscopic variables. Fig 3E  shows all three key observables and shows that the DME method is in an excellent agreement with the model.
In the previous example we demonstrated that the method works well in the most extreme case, when the forces 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 may suggest that the method will perform even better under slow environmental changes, which in reality may be 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 environmental changes can be observed on various timescales, for example the warming of the oceans [34] or the yearly cycle of seasons [35], resulting in different migration rates throughout the year due to varying resource abundance.
In Fig 4 we compare periodic changes for three scenarios: an abrupt, a slow and a fast but continuous change. We 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. While the faster environmental oscillations result in a smaller range of average population sizes, the range of effective forces, on the contrary, increases. In the extreme case of very fast environmental changes (e.g. oscillations with large frequency, or environmental temporal noise of a fixed variation) one expects that the average population size will stay almost stationary as the convergence to equilibrium is much slower than the timescale of environmental fluctuations, while the rapid fluctuations will impact the effective forces, similarly to our first scenario in Fig 4A. We found that the solution of DME is in a good agreement with that of the FPE, with the temporal dynamics in the DME and the moments in the FPE equation being indistinguishable by eye (white or black lines). We quantified the error of the DME approximation by the relative entropy between the exact population density, which solves the FPE equation, and the DME approximation, which is computed from the dynamics of the effective forces. The magnitude of this error is larger in the scenario A where the environmental parameters change abruptly, as compared to their smooth changes in cases B-C (note the diferences in the range between the panels). Surprisingly, even in the first scenario where the environmental parameters change abruptly, and the dynamics have time to adapt to the changed environment, the relative entropy does not decay monotonically during the adaptation period. We observe that after a brief initial exponential drop (showed on a logarithmic scale in S1 Fig), the relative entropy follows a slower non-monotonous pattern. This is likely caused by the dynamics of the effective forces, each of which shows non-monotonous convergence to the new equilibrium at a different time scale. This non-monotonicity of the relative entropy, which we observe in all studied scenarios, makes the problem difficult to study analytically.

Discussion
We presented an application of the dynamic maximum entropy method, which helps reduce complexity of the stochastic process by linking microscopic quantities to macroscopic observables.
We first studied the Ornstein-Uhlenbeck process with time-dependent parameters. We used the DME method-instead of following the full stochastic dynamics, we followed just the key observables (in this case the first two moments), which change deterministically. The observables and the forces acting on them were identified from the potential form of the stationary solution (in this case Gaussian). The dynamic problem was solved by the DME method, which uses the stationary ansatz with time-dependent parameters to best approximate the dynamics of observables. We derived a two-dimensional dynamical system of the effective forces that characterize the solution of the OU dynamics. We showed that despite the intricacies of the DME approximation, the DME dynamics coincide with an exact solution of the OU process. This is because the dynamical equations for the first two moments are closed and thus the solution is Gaussian at all times. However, even though the OU dynamics are linear, the effective forces solve a nonlinear system of ordinary differential equations.
The focus of our work is on the stochastic island model, represented by nonlinear population dynamics based on a logistic equation supplemented by migration from other islands.
The key parameters of the problem, the intrinsic growth rate, the carrying capacity, and the migration rate, are in general functions of time, reflecting temporal environmental changes, which take the problem out-of-equilibrium. Nonlinearity of the process results in the dynamics of moments hn k i which are not closed for any k. Therefore we used DME to derive a 3-dimensional nonlinear dynamical system for the effective forces using DME. The associated observables are no longer just the first three moments but include hni, h−n 2 /2i and hlogni. Unlike for the OU model, the DME approximation of the island model with migration is no longer exact. The system is not fully explicit but contains terms using hypergeometric functions. Nevertheless, it can be solved for dynamic environmental 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 changes of the environmental forces the observables in the DME approximation are still extremely accurate at all times and thus the effective forces serve as a proxy for the dynamics. When the environmental forces settle to constant values, the effective forces also converge to these values. The DME serves as a change of optics where we represent the non-equilibrium dynamics using a series of equilibria parametrized by dynamical effective forces. The strength of the DME method in our view is not in speeding up the numerical method for solving the problem but in better understanding the underlying dynamics in an appropriate low-dimensional space.
Although the FP equations give a complete description of how the whole probability distribution evolves through time, it is not feasible to solve them numerically for more than two or three variables, and analytic results become intractible, once one moves away from linear models. The FP equations do not bring much actual understanding, without the help of heuristic approximations such as ME. Although it is possible to use the quasi-stationary approximation without the connection to ME, the method provides us with useful information. By formulating the suitable ME problem for the stationary distribution we learned which macroscopic quantities are important for the low-dimensional projection of the dynamics. More specifically, DME identifies pairs of complementary variables-for example, in the population dynamic example, (n, −n 2 , log(n)) correspond to growth rate, strength of density regulation, and migration. While the observables enter the variational problem via the constraints, the forces are the corresponding Lagrange multipliers in ME. This is an intriguing extension to more general nonlinear stochastic processes of an idea familiar from statistical thermodynamics.
In addition, the DME method can be placed into the arsenal of methods for non-equilibrium dynamics, as we have shown in Fig 1. It is built from the stationary ME method but unlike the MC method, which is essentially a ME method applied on temporal trajectories, it uses the FPE to establish relationships between the time points.
One of the most striking properties of the DME method is its accuracy on the macroscopic level. This is surprising because the quasi-stationary assumption suggests validity of the approach only when the applied forces change slowly (i.e., adiabatically). However, even for fast-changing forces the approximation stays very accurate. We have shown numerically that the relative entropy between the exact and approximate solutions does not decay monotonically even in the simplest scenario of an abrupt change in the ecological forces. This along with the unusual form of the DME approximation makes analytical study of the accuracy of the method a difficult mathematical problem which remains open to date, despite insight provided in [21].
The remarkable accuracy of the DME approximation, even when conditions change abruptly, remains mysterious to us. We can draw an analogy with path integration, where the approximation that fluctuations around the most likely path are Gaussian is exactly correct for all cases that can be solved explicitly, and is accurate even with strong nonlinearity [36]. The breakdown of DME when m < 1 2 is more understandable: then, the concentration of probability near the small population boundary in effect constrains one of the free variables, so that the approximation cannot fit as well, based on the remaining variables. This problem can be resolved by considering boundary layer in the probability distribution, thus effectively increasing the number of free variables, as investigated in detail in [16].
This work outlines the first step towards studying more complex questions where the complexity of the 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 interaction has been studied in [32] but only in the stationary case. Although the existence of an explicit stationary distribution 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 on the structure of the problem, in which the stochastic dynamics are described by the FPE and the stationary solution is explicit. This includes a wide range of problems accross disciplines, for example stochastic coagulation-fragmentation dynamics when the rates satisfy a detailed balance condition (existence of an explicit stationary distribution for this problem was shown in [37]).