Exact dynamics of stochastic linear delayed systems: Application to spatiotemporal coordination of comoving agents

Delayed dynamics result from finite transmission speeds of a signal in the form of energy, mass, or information. In stochastic systems the resulting lagged dynamics challenge our understanding due to the rich behavioral repertoire encompassing monotonic, oscillatory, and unstable evolution. Despite the vast literature, quantifying this rich behavior is limited by a lack of explicit analytic studies of high-dimensional stochastic delay systems. Here we fill this gap for systems governed by a linear Langevin equation of any number of delays and spatial dimensions with additive Gaussian noise. By exploiting Laplace transforms we are able to derive an exact time-dependent analytic solution of the Langevin equation. By using characteristic functionals we are able to construct the full time dependence of the multivariate probability distribution of the stochastic process as a function of the delayed and nondelayed random variables. As an application we consider interactions in animal collective movement that go beyond the traditional assumption of instantaneous alignment. We propose models for coordinated maneuvers of comoving agents applicable to recent empirical findings in pigeons and bats whereby individuals copy the heading of their neighbors with some delay. We highlight possible strategies that individual pairs may adopt to reduce the variance in their velocity difference and/or in their spatial separation. We also show that a minimum in the variance of the spatial separation at long times can be achieved with certain ratios of measurement to reaction delay.

Figure 1

Graphical representation of the various steps involved in solving a 1D linear stochastic delay differential equations with $N$ delays, with corresponding rates ${\alpha }_i$ and noise strength $\sigma$. First, the propagator $\Lambda \left(t\right)$ of the associated deterministic differential equation of the system (top box) is found by solving the transformed equation in the Laplace domain, i.e., $\tilde{\Lambda }\left(\epsilon \right)$, and then inverting it back to the time domain. The propagator is found to have two functional representations: an infinite sum of exponentials (when all poles in Laplace space are simple) or a sum of delay shifted polynomials. With the time dependence of the propagator it is then possible to build a multivariate probability distribution for the dynamics of the system as well as to obtain the mean-square displacement (MSD) at finite time, via integration of ${\Lambda }^2\left(t\right)$, and at infinite time through a Fourier inversion ${\mathcal{F}}^{-1}$ of the modulus square of the Fourier transform of the propagator $\widehat{\Lambda }\left(\omega \right)$. For simplicity we have presented here only the 1D case, but the same sequential steps can be followed when the system possesses more than one degree of freedom.

Figure 2

Temporal dynamics of the two-delay propagator $\Lambda \left(t\right)$, governed by the deterministic delay differential equation $\dot{\Lambda }\left(t\right)=-\Lambda (t-{\tau }_1)-\Lambda (t-{\tau }_2)$ with the dimensionless delays ${\tau }_1=1/3$, ${\tau }_2=1/2$, plotted from Eq. (5). The number of required summation terms is also plotted as a function of time through the steplike solid blue line. Each jump of the blue curve indicates when the analytic solution changes functional form, corresponding to when updates in the rate of change of the differential equation occurs. The magnification window shows how the summation terms contribute at short times, with $k_1$ and $k_2$ marked explicitly. The corresponding ordered update times are at $t=1/3$, $t=1/2$, $t=2/3$, $t=5/6$, $t=1$, and $t=7/6$. At each time $t$ the order of the polynomial function in Eq. (5) can be evinced by the value of $k_1+k_2$, e.g., at times $t=5/6$, 1, and 7/6, it is, respectively, of orders 2, 5, and 3.

Figure 3

Dynamics and variance of the one delay equation $\dot{x}\left(t\right)=-\alpha x(t-\tau )+\sigma \xi \left(t\right)$ with $\varphi \left(t\right)=x\left(0\right)$ for $-\tau \le t<0$. Panel (a) displays a contour plot of the long-time MSD. The plot shows the variation with the delay $\tau$ and rate $\alpha$. The parameter region for instability $\alpha \tau \ge \pi /2$ and oscillations $e^{-1}\le \alpha \tau <\pi /2$ are marked along with the minimum of the variance, given by $\alpha \tau =cos\left(\alpha \tau \right)$. In panel (b) sample propagator dynamics are plotted for a spectrum of $\alpha =\tau$ parameter values and the corresponding trajectories of $x\left(t\right)$ are plotted in panel (c) for the initial condition $x\left(0\right)=1$. The mean-squared displacement corresponding to these trajectories is plotted in panel (d). The difference in the dynamics between $\langle x\left(t\right)\rangle$ and the propagator in this example is simply a shift in time with respect to the delay because $x\left(t\right)=-\alpha x(t-\tau )$ with $x\left(t\right)=1$ for $-\tau \le t\le 0$, whereas $\Lambda \left(t\right)=-\alpha \Lambda (t-\tau )$ with $\Lambda \left(0\right)=1$ and $\Lambda \left(t\right)=0$ for $-\tau \le t<0$.

Figure 4

Stationary MSD for the two-delay system governed by $\dot{x}\left(t\right)=-\alpha [x(t-{\tau }_1)+x(t-{\tau }_2)]+\sigma \xi \left(t\right)$. Panel (a) shows the variation in MSD, for fixed $\alpha =1$, as a function of the smallest delay ${\tau }_1$ and the delay difference ${\tau }_2-{\tau }_1$, found by numerically integrating Eq. (9). The MSD for the nondelayed case ${\sigma }^2/4\alpha$ (Ornstein-Uhlenbeck) has been subtracted everywhere from the delayed MSD to highlight the difference between these cases. The stability boundary [83] is also shown and the MSD is seen to diverge approaching this boundary. Panel (b) shows the variation in MSD as a function of the coefficient $\alpha$ for four specific cases and a spectrum of delay values. The markers show the results of numerically integrating and ensemble averaging the Langevin equation while the solid lines have been found analytically using Eq. (9).

Figure 5

Timeline showing the delayed nature of the reaction of individual $i$, at time $t$, to the action of individual $j$, at time $t-{\tau }_r-{\tau }_m$. The delay between action and reaction occurs because it takes a finite amount of time, ${\tau }_m$, for individual $i$ to measure the changes in individuals $j$'s properties and a finite amount of time for individuals $i$ and $j$ to decide and perform a movement change. Individual $i$ will thus produce a movement change at time $t$ with a delay ${\tau }_r$ and a delay ${\tau }_r+{\tau }_m$ with respect to the movement changes of individual $j$.

Figure 6

Representation of the dynamics of the velocity matching interaction. In panel (a) we show trajectories from the model in Eq. (28) with measurement and response delays given in the legend. In the inset the velocity separations for the deterministic component of the trajectories in the $x$ and $y$ directions are also plotted. In panel (b) we present a scatter plot of the dynamics of the velocity separation in real time and delay time for the case ${\tau }_m=0.4$, ${\tau }_r=0$ presented in panel (a). In the left and bottom adjoining plots we show the corresponding marginal probability densities at four specific times, $t/\tau =1/4$, $1/2$, $3/4$, and 1, with the $x$ and $y$ components of the velocity plotted as circles and squares, respectively. The $\sigma$ ellipses for the stationary multivariate probability in the main figure have been plotted from the associated Fourier expressions in (16). The time-dependent probabilities for the left and bottom adjoining panels are calculated by marginalizing, respectively, the delay and nondelay component from the full multivariate probability expression derived from Eq. (26).

Figure 7

Dynamics of higher-order interaction coordination. Three sample alignment trajectories from simulation of Eq. (34) are plotted in panel (a) with the parameters given in the legend. In the inset figures the corresponding speed and distance separations are displayed. The dark blue curve does not converge to the expected value $0.1\sqrt{2}\simeq 0.14$ since it is purely a velocity matching interaction, i.e., $\beta =0$ in Eq. (34). In panel (b) we draw a scatter plot of the spatial versus velocity separation for the middle trajectory in (a), with adjoining plots showing the respective marginal probability densities at three specific times, $t/{\tau }_x=25/8$, $30/8$, and $35/8$. The $x$ and $y$ directions are plotted as circles and squares, respectively. The $\sigma$ ellipses for the stationary multivariate probability have been plotted as well as the time-dependent probabilities for the marginal distributions.

Figure 8

Contour level maps of the velocity separation variance in (a), and the spatial separation variance in (b), as a function of the coupling ratio, $\rho =\beta /\alpha$, and delay ratio, $\eta ={\tau }_x/{\tau }_v$. The plots are found by numerically integrating the equations for the stationary covariance up to a divergence limit used to mark the instability boundary. The velocity variance is seen to increase monotonically towards the instability boundary, whereas the spatial variance contains a minimum with respect to the coupling ratio, shown by the white line.

Figure 9

Eigenvalues, $\lambda$, of the system $\dot{x}\left(t\right)=-\alpha x(t-\tau )$, given by (A8), calculated using the first 40 branches of the Lambert function for a spectrum of parameter values $\alpha \tau$. The values where oscillations and instability are induced have been highlighted in black.

Figure 10

Stationary variance for the mutual velocity matching interaction (29) as a function of the coupling strength $\gamma$ and the delay ratio $\rho ={\tau }_r/({\tau }_r+{\tau }_m)$ with fixed ${\tau }_r=1/2$. The variance diverges approaching the instability boundary in the top right, whereas a minimum, displayed by the white line, appears as a function of $\rho$. The white circles show the location of this minima as approximated by Hunt et al. [83], which we have reported for completeness in Eqs. (D2) and (D3).