Deterministic particle flows for constraining stochastic nonlinear systems

Devising optimal interventions for constraining stochastic systems is a challenging endeavor that has to confront the interplay between randomness and dynamical nonlinearity. Existing intervention methods that employ stochastic path sampling scale poorly with increasing system dimension and are slow to converge. Here we propose a generally applicable and practically feasible methodology that computes the optimal interventions in a noniterative scheme. We formulate the optimal dynamical adjustments in terms of deterministically sampled probability flows approximated by an interacting particle system. Applied to several biologically inspired models, we demonstrate that our method provides the necessary optimal controls in settings with terminal, transient, or generalized collective state constraints and arbitrary system dynamics.

Figure 1

Schematic of forward and time-reversed probability flows for deriving state- and time-dependent dynamical interventions $u(x,t)$. We initially sample the flow ${\varrho }_t\left(x\right)$ for the time interval $[0,T]$. By employing the logarithmic gradient (score) of ${\varrho }_t\left(x\right)$, we evolve the time-reversed constrained probability flow ${\tilde{q}}_{\tau }\left(x\right)$. The optimal state- and time-dependent dynamical interventions $u(x,t)$ result from the difference of the logarithmic gradients of the two probability flows.

Figure 2

Deterministic particle flow (DPF) control provides optimal interventions to drive the controlled system to target state (grey cross). (a) A controlled trajectory starting from state ${\mathbf{x}}_0=(-1,1)$ reaches the target ${\mathbf{x}}^{*}=(1,1)$ at time $T=0.7$ (blue-yellow line), while an uncontrolled trajectory remains in the vicinity of initial condition ${\mathbf{x}}_0$ during the same time interval (orange line). (b) Mean and standard deviation of the marginal densities of 1000 controlled trajectories employing interventions computed with our framework DPF (purple) and with PICE (grey). Orange line indicates mean of 1000 uncontrolled trajectories, while shaded area captures the associated standard deviation. For estimating the controls, we employed $N=400$ particles for DPF, and $N_{\mathrm{pice}}=500$ for PICE. (c.) Comparison of (upper) (logarithmic) control energy ${\parallel u(\mathbf{x},t)\parallel }_2^2$, and (lower) deviation from the terminal point $\parallel {\mathbf{X}}_T-{\mathbf{x}}^{*}{\parallel }^2$ for each controlled trajectory with interventions computed according to DPF (magenta) and PICE (grey). (d) (Logarithmic) control energy (upper) improves moderately, and terminal error (lower) remains constant for increasing particle number $N$. The number of inducing points in the logarithmic gradient estimation conveys negligible difference in control energy and terminal error (inducing point number magenta: $M=50$, green: $M=100$). Grey line indicates the performance of PICE in the same setting. [(e)–(h)] Same as (a)–(d) with additional path constraint $U((x,y),t)=\beta {(y-1)}^2$, with $\beta ={10}^3$. Average iteration number required for PICE was 250 in the non path-constrained setting, and 300 for the path constrained one.

Figure 3

Equivalence of dynamical interventions delivered by deterministic particle flow control (DPF), guided deterministic particle flow control (gDPF), and PICE, for typical (top) and atypical (bottom) terminal conditions. (Top) (a) Controlled trajectories simulated by employing dynamical interventions delivered by deterministic particle flow control (magenta), guided deterministic particle flow control (orange), and path integral cross entropy method (grey). All trajectories started from initial point $x_0=-1$ (left silver circle) and reached the target $x^{*}=0$ at time $T=1$. (b) Transient mean ${\mu }_t$ and standard deviation ${\sigma }_t$ over 1000 independent trajectories controlled with each framework. (c) Total control energy ${\parallel u(x,t)\parallel }^2$, and (d) deviation from terminal point for 1000 independent controlled trajectories with interventions computed according to each framework. Grey horizontal lines in (c) and (d) denote mean values over all realizations. The proposed methods result in slightly more expensive control costs, but are more consistent in precisely reaching the terminal state. (Bottom) Same as upper row for target $x^{*}=1$ only for gDPF and PICE. Here DPF is not applicable since the forward probability flow does not reach the atypical target point $x^{*}=1$.

Figure 4

Optimal interventions effectively drive the nonconservative system to the target state for different noise amplitudes. (a) Individual trajectory controlled by DPF (green-yellow) starting from stable state $x_0=(1.996,0)$ successfully reaches the target $x^{*}=(1,1)$ at $T=0.5$, while the uncontrolled trajectory (grey) fails to leave the basin of attraction of $x_0=(1.996,0)$. (b) Agreement between controlled densities of 1000 independent controlled trajectories driven by DPF (magenta) and PICE (grey) for noise amplitude $\sigma =1.5$. (c) Dissipated control energy and (d) terminal error for increasing noise strength $\sigma$ for the two control frameworks. For increasing noise DPF delivers more efficient, but moderately less precise controls than PICE. Further parameters: particle number: $N=600$ inducing point number: $M=20$, discretization time step $dt={10}^{-3}$.

Figure 5

Deterministic particle flow (DPF) control steers an evolutionary process of molecular phenotypes to a non (evolutionary) optimum state. [(a) and (b)] Transient mean of marginal densities of 1000 controlled trajectories, ${\mu }_x^{q_t}$ and ${\mu }_y^{q_t}$, of controlled trajectories of an evolutionary process with covarying phenotypes with correlation (a) ${\varrho }_{xy}=0.25$, and (b) ${\varrho }_{xy}=0.75$ starting from phenotypic state ${\mathbf{x}}_0=(0,0)$ and reaching the target ${\mathbf{x}}^{*}=(0.5,0)$ (green cross) at time $T=1.0$. The employed interventions were computed with the correct model (solid lines) or with a misspecified model that assumed zero cross-phenotype correlations ${\varrho }_{xy}=0.$ (dashed line) imposing: (i) only terminal constraints (purple lines) and (ii) both terminal and path constraints (light blue lines lines). The path constraint $U((x,y),t)=\beta y^2$, with $\beta ={10}^3$ was set to limit the fluctuations around the $y$ axis. The white line and grey shaded region indicate the mean and standard deviation of 1000 uncontrolled trajectories evolving under natural selection from the same initial state ${\mathbf{x}}_0$. [(c) and (d)] Comparison of (c) dissipated control energy ${\parallel u(x,t)\parallel }_2^2$, and (d) path deviations from target along the second ($y$) phenotype axis $\int ({\left(X_t\right)}_y-y^{*})dt$ when steering the evolutionary processes with (top) anticorrelated and (bottom) correlated phenotypes with increasing (in absolute value) correlation coefficient ${\varrho }_{xy}$ from left to right, with ${\varrho }_{xy}=\pm \{0.25,0.50,0.75\}$. The employed interventions followed the same model and constraint assumptions as in (a) and (b) with “cor” and “mis” denoting results from computed controls employing the correct and the misspecified dynamical model, respectively. Further parameters: particle number: $N=1000$, inducing point number: $M=50$, discretization time step: $dt={10}^{-3}$.

Figure 6

Synchronization control of two coupled Kuramoto phase oscillators. (a) Evolution of phases (${\theta }_i$) of two controlled (purple) and two uncontrolled (green) Kuramoto oscillators mutually coupled with coupling $J=1.2$ and noise amplitude $\sigma =1.0$. (b) Evolution of Kuramoto phase-coherence order parameter $R$ for the controlled (purple) oscillators indicates a fast transition to synchrony ($R=1$), while the identical uncontrolled oscillators become progressively incoherent, indicated by a strongly fluctuating order parameter (green). The orange line denotes the expected long time average value of the order parameter for noninteracting oscillators considering finite size scaling effects. The grey line marks the level of $R=1$ indicating a completely synchronous state. (c) Control input provided to each oscillator. Further parameters: particle number: $N=2000$, inducing point number: $M=80$, natural frequencies: ${\omega }_i=0$, initial condition: ${\theta }^{\left(i\right)}\sim \mathcal{N}(3,0.5^2)$, and $T=1.5$ (time units).

Figure 7

Synchronization control of a network of two coupled Kuramoto phase oscillators for different coupling strengths. (a) Time averaged phase-coherence order parameter ($R$) of controlled (purple) and uncontrolled (green) networks under different noise conditions [$\sigma =0.5$ (triangles) and 1.0 (circles)]. The proposed method (DPF) effectively synchronizes the controlled oscillators already for vanishing coupling strength between them. (b) Evolution of Kuramoto order parameter $R$ for networks with coupling $J=2.0$ and two noise conditions [$\sigma =0.5$ (top) and 1.0 (bottom)] for controlled (purple) and uncontrolled (green) oscillators. The control induces fast transition to synchrony ($R=1$), while the identical uncontrolled oscillators either synchronize slower (for low noise), or become only partially synchronized (for strong noise). Individual lines indicate evolution of the order parameter in 20 realizations of the network starting from same initial conditions and from a single computation of the required controls for each setting (where relevant). Dotted black lines denote the mean over the 20 realizations. Further parameters: particle number $N=2000$, inducing point number: $M=80$, natural frequencies: ${\omega }_i=0$, initial condition: ${\theta }^{\left(i\right)}\sim \mathcal{N}(3,0.5^2)$ and $T=1.5$ (time units).

Figure 8

Onset of synchronization and percentage of time in synchronized state after synchrony onset reveal the effectiveness of deterministic particle flow control to induce robust synchronization on a network of $K=2$ oscillators. (a) Onset of synchronization for controlled (purple) and uncontrolled (green) networks quantified as the first time $t^{\mathrm{syn}}$ the phase-coherence order parameter exceeds $R\ge 0.99$ and remains over that value for duration ${\tau }^s=20\times dt=0.02$ time units. Grey annotations denote the percentage of the examined networks that reached the synchronous state for duration ${\tau }^s$. Absence of annotation indicates that all examined networks reached synchrony. (b) Percentage of simulation time the networks spontaneously spent in synchronized state ($R\ge 0.99$) after synchrony onset $t^{\mathrm{syn}}$. Both figures consider two different noise conditions [$\sigma =0.5$ (triangles) and 1.0 (circles)]. Further parameters: particle number $N=2000$, inducing point number: $M=80$ and $T=1.5$ (time units). For each noise condition dots denote average over 3 control computations with different initial conditions with 20 controlled trajectories for each (60 total controlled trajectories for each point).

Figure 9

Synchronization control of a finite-size network of $K=6$ interacting Kuramoto phase oscillators. (a) Evolution of phases ${\theta }^{\left(i\right)}$ of a controlled network and (b) an identical uncontrolled network. The phases of the oscillators quickly synchronize when controlled by DPF and remain synchronized throughout the entire simulation interval $[0,T=0.5]$. In the absence of control the phases of the oscillators become increasingly incoherent. (c) Evolution of Kuramoto order parameter capturing phase coherence $R$ for the controlled (purple) oscillator network indicates a rapid transition to complete synchrony ($R=1$) for all 20 realizations (individual purple lines). The order parameter of identical uncontrolled networks fluctuates strongly indicating partial incoherence (green). The orange line denotes the expected long time average value of the order parameter if the oscillators were noninteracting considering finite size scaling effects. For visual clarity, the grey line marks the level of $R=1$ indicating a completely synchronous state. (d) Control energy spent on all $K=6$ oscillators for a single control realization. Further parameters: coupling strength $J=1.0$, noise amplitude $\sigma =1$, particle number $N=3000$, inducing point number $M=300$, initial condition: ${\theta }^{\left(i\right)}\sim \mathcal{N}(3,0.5^2)$, natural frequencies: ${\omega }_i\sim \mathcal{N}(0,1)$ discretization time step $dt={10}^{-3}$ and $T=0.5$ (time units).

Figure 10

Synchronization control of a network of six heterogeneous Kuramoto phase oscillators for different coupling strengths. (a.) Time averaged phase-coherence order parameter ($R$) of controlled (purple) and uncontrolled (green) networks under different noise conditions [$\sigma =1.0$ (triangles) and 1.5 (circles)]. The proposed method (DPF) effectively synchronizes the controlled oscillators already for weak coupling. (b.) Evolution of Kuramoto order parameter $R$ for networks with coupling $J=3.0$ and two noise conditions [$\sigma =1.0$ (top) and 1.5 (bottom)] for controlled (purple) and uncontrolled (green) oscillators. The control induces fast transition to synchrony ($R=1$), while the identical uncontrolled oscillators either synchronize slower (for low noise), or become only partially synchronized (for strong noise). Individual lines indicate evolution of the order parameter in 20 realizations of the network starting from same initial conditions and from a single computation of the required controls for each setting (where relevant). Dashed black lines denote the mean over the 20 realizations. Further parameters: particle number $N=3000$, inducing point number: $M=300$, initial condition: ${\theta }^{\left(i\right)}\sim \mathcal{N}(3,0.5^2)$, natural frequencies: ${\omega }_i\sim \mathcal{N}(0,1)$ and $T=0.5$ (time units). Each point in (a) denotes the mean over 20 realizations of a single control computation.