Exploring the Ant Mill: Numerical and Analytical Investigations of Mixed Memory-Reinforcement Systems

TL;DR

This paper models the complex motion of ant swarms, especially the death spiral phenomenon, using a novel continuous random walk framework that combines memory and reinforcement effects, supported by analytical and numerical stability analysis.

Contribution

It introduces a new continuous PDE-based model integrating memory and reinforcement in random walks, capturing the death spiral phenomenon in ant swarms.

Findings

Derived an axi-symmetric, steady-state solution resembling the death spiral.

Proved the stability of the steady-state solution numerically.

Demonstrated the model's ability to replicate complex biological motion patterns.

Abstract

Under certain circumstances, a swarm of a species of trail-laying ants known as army ants can become caught in a doomed revolving motion known as the death spiral, in which each ant follows the one in front of it in a never-ending loop until they all drop dead from exhaustion. This phenomenon, as well as the ordinary motions of many ant species and certain slime molds, can be modeled using reinforced random walks and random walks with memory. In a reinforced random walk, the path taken by a moving particle is influenced by the previous paths taken by other particles. In a random walk with memory, a particle is more likely to continue along its line of motion than change its direction. Both memory and reinforcement have been studied independently in random walks with interesting results. However, real biological motion is a result of a combination of both memory and reinforcement. In this paper, we construct a continuous random walk model based on diffusion-advection partial differential equations that combine memory and reinforcement. We find an axi-symmetric, time-independent solution to the equations that resembles the death spiral. Finally, we prove numerically that the obtained steady-state solution is stable.

Equations

$$\frac{d\rho }{dt}=\frac{\partial \rho }{\partial t}+\vec{v}\cdot {\nabla }_x\rho =-\alpha \rho +\alpha \int T(\vec{v},{\vec{v}}^{\prime })\rho (\vec{x},{\vec{v}}^{\prime },t)d{v}^{\prime },$$

$$\frac{d\rho }{dt}=\frac{\partial \rho }{\partial t}+\vec{v}\cdot {\nabla }_x\rho =-\alpha \rho +\alpha {\int }_{-\pi }^{\pi }T(\theta ,{\theta }^{\prime })\rho (\vec{x},{\theta }^{\prime },t)d{\theta }^{\prime }.$$

$$\frac{\partial g}{\partial t}=\gamma {\int }_{-\pi }^{\pi }\rho (x,{\theta }^{\prime },t)d{\theta }^{\prime }-\beta g.$$

$$T(\theta ,{\theta }_g)=\frac{1}{4{\pi }^2}(Jcos(\theta -{\theta }_g)+1),$$

$$\frac{\partial \rho }{\partial t}+v(\cos \theta ,\sin \theta )\cdot {\nabla }_x\rho =-\alpha \rho +\frac{\alpha }{4{\pi }^2}(Jcos(\theta -{\theta }_g)+1){\int }_{-\pi }^{\pi }\rho (x,{\theta }^{\prime },t)d{\theta }^{\prime }.$$

$$(v\cdot 2\pi ik\cos \theta +\alpha )P(ik,\theta )-\frac{\alpha }{4{\pi }^2}(Jcos\theta +1){\int }_{-\pi }^{\pi }P(ik,{\theta }^{\prime })d{\theta }^{\prime }=0.$$

$$\frac{\partial \rho }{\partial t}=D\frac{\partial }{\partial x}\left(\frac{\partial \rho }{\partial x}-\rho \frac{\beta }{\alpha +\beta g}\frac{\partial g}{\partial x}\right)$$

$$\frac{\partial g}{\partial t}=\lambda \rho -g$$

$$\frac{d\rho }{dt}=\frac{\partial \rho }{\partial t}+\vec{v}\cdot \nabla \rho =\nabla \cdot (\nabla \rho -\rho \frac{\beta }{\alpha +\beta g}\nabla g).$$

$$\frac{d\vec{v}}{dt}=\frac{\partial \vec{v}}{\partial t}+\vec{v}\cdot \nabla \vec{v}=b\nabla g$$

$$\frac{d\rho }{dt}=\frac{\partial \rho }{\partial t}+\vec{v}\cdot \nabla \rho =\nabla \cdot (\nabla \rho -\rho \frac{\beta }{\alpha +\beta g}\nabla g).$$

$$\frac{\partial g}{\partial t}=\lambda \rho -g$$

$$\frac{d\vec{v}}{dt}=\frac{\partial \vec{v}}{\partial t}+\vec{v}\cdot \nabla \vec{v}=b\nabla g.$$

$$\frac{d\rho }{dt}=\frac{\partial }{\partial r}\left(r\frac{\partial \rho }{\partial r}-r\rho \frac{\beta }{\alpha +\beta g}\frac{\partial g}{\partial r}\right)$$

$$\frac{\partial g}{\partial t}=\lambda \rho -g$$

$$\frac{d\vec{v}}{dt}=\frac{\partial \vec{v}}{\partial t}-\frac{v_{\theta }^2}{r}=b\frac{\partial g}{\partial r}$$

$$0=\frac{\partial }{\partial r}\left(r\frac{\partial \rho }{\partial r}-r\rho \frac{\beta }{\alpha +\beta g}\frac{\partial g}{\partial r}\right)$$

$$0=\lambda \rho -g$$

$$-\frac{v_{\theta }^2}{r}=b\frac{\partial g}{\partial r}.$$

$$g=\lambda \rho$$

$$0=\frac{\partial }{\partial r}\left(r\frac{\partial \rho }{\partial r}-r\lambda \rho \frac{\beta }{\alpha +\beta \lambda \rho }\frac{\partial \rho }{\partial r}\right).$$

$$\int 0dr=\int \frac{\partial }{\partial r}\left(r\frac{\partial \rho }{\partial r}-r\lambda \rho \frac{\beta }{\alpha +\beta \lambda \rho }\frac{\partial \rho }{\partial r}\right)dr$$

$$C_1=r\frac{\partial \rho }{\partial r}-r\lambda \rho \frac{\beta }{\alpha +\beta \lambda \rho }\frac{\partial \rho }{\partial r}$$

$$C_1(\alpha +\beta \lambda \rho )=r\frac{\partial \rho }{\partial r}(\alpha +\beta \lambda \rho )-r\lambda \rho \beta \frac{\partial \rho }{\partial r}$$

$$C_1(\alpha +\beta \lambda \rho )=r\alpha \frac{\partial \rho }{\partial r}$$

$$r\frac{\partial \rho }{\partial r}+C_1\frac{\beta }{\alpha }\lambda \rho =-C_1$$

$$r\frac{\partial \rho }{\partial r}=-C_1-C_1\frac{\beta }{\alpha }\lambda \rho$$

$$\frac{1}{-C_1-C_1\frac{\beta }{\alpha }\lambda \rho }\partial \rho =\frac{1}{r}\partial r$$

$$\int \frac{1}{-C_1-C_1\frac{\beta }{\alpha }\lambda \rho }\partial \rho =\int \frac{1}{r}\partial r$$

$$-\frac{1}{C_1}\int \frac{1}{1+\frac{\beta }{\alpha }\lambda \rho }\partial \rho =\int \frac{1}{r}\partial r$$

Full text

1 Introduction

The process by which living organisms respond to chemical stimuli in their environments by moving towards or away from them is called chemotaxis. Movement towards a higher concentration of chemical is known as positive taxis (taxis in Greek means to arrange), while movement towards a lower concentration of chemical is known as negative taxis (Friedman and Tello, 2002). Chemotaxis has been observed in several species of bacteria, including the well-studied soil species known as the myxobacteria; slime molds; and ant species such as army ants (Othmer and Stevens, 1997). A chemotaxis process also occurs in the growth of new blood vessels, known as angiogenesis, that occurs in the expansion of a tumor. In tumor-induced angiogenesis, chemical agents secreted by a tumor attract neighboring endothelial cells, which form the surface of blood vessels that then sprout towards the tumor to provide nourishment (Friedman and Tello, 2002).

The mathematics that underlies each of these chemotaxis scenarios is based on an object known as the random walk. A random walk is any path that consists of a sequence of random steps. Examples of random walks include the path traveled by a molecule in a liquid or gas and the path of a foraging animal. The class of random walks that have been used in the past to model chemotaxis processes are called reinforced random walks (RRWs). In this type of random walk, the walker can modify (reinforce) the favorability of traveling along a particular path, for example by dropping an attracting chemical substance. A variety of different reinforced random walk systems based on differing reinforcement mechanisms have modeled many observed chemotaxis situations with success (Pemantle, 2006).

Although reinforced random walk models have been able to replicate several different chemotaxis phenomena, there remain limitations to these systems (Codling, Plank, and Benhamou, 2008). In particular, reinforced random walks do not take into account another important property of animal and cell motion: the tendency of these bodies to continue moving in the same direction. This property of inertia is explicitly considered in a class of random walks called correlated random walks (CRWs). In these walks, there is a correlation between successive step orientations, which produces a local directional bias that causes each step to tend to point in the same direction as the previous one. The influence of the initial direction of motion gradually diminishes over time.

The motion of real biological systems is complex and is based on both correlation, which is also termed memory, and reinforcement. The fact that these two factors influencing motion have only been studied independently may explain why there continue to be observed phenomena that have not been replicated by past models. One such phenomenon is the army ant death spiral. In this example, we suspect that there are two basic directional forces influencing ant motion: $(1)$ a chemical-gradient-induced centripetal force and $(2)$ a tangential force related to inertia. These two influences are shown for a single ant caught in the spiral in Figure 2. The concentration of chemical is highest in a ring around the center of the spiral and fades out radially, creating a gradient that points toward the center for every point that the ants are located. The “eye” of the spiral corresponds to the region where no ants are located in Figure $1$. The combined effect of this centripetal force and the tangential pull of inertia produces the observed circular motion of every ant in the configuration.

If reinforcement based on the trails laid by the ants was the only factor influencing their motion, it is unlikely that they would maintain circular motion. Instead, they would all move closer to the center of the spiral, where the chemical concentration was higher. Similarly, if the chemical gradient had a negligible effect, and inertia was the only bias propelling the ants, they would also not have circular motion. We seek a mathematical description that finds a balance between these two forces so as to replicate the army ant spiral that we observe in real life.

In this paper, we develop mixed memory-reinforcement models with the aim of replicating the army ant death spiral. Moreover, we aim to prove the stability of the ant death spiral mathematically based on the developed models. We begin our exploration of mixed memory and reinforcement in random walks with intuitive models based on integral equations (Section $2$). Due to the determined limitations of these systems of equations, we then consider random walk models based on diffusion-advection equations. We construct a system of diffusion-advection partial differential equations that has a time-independent solution resembling the death spiral, and then numerically prove the solution is stable in time (Section $3$).

2 Integral equation model

We begin our quest for a mathematical model of the ant mill by developing an intuitive system that describes all properties of ant motion as simply as possible. In this model, we let the position $\vec{x}$ of an ant in a $2$-dimensional plane be independent of the ant’s velocity, $\vec{v}$. For every $\vec{x}$ in position space, every $\vec{v}$ in velocity space, and every time $t>0$, the density $\rho(\vec{x},\vec{v},t)$ is the density of ants at that point in the position-velocity phase space at time $t$. The system thus keeps track of position and velocity information for every ant in its domain over time. We can simplify this model further by placing the additional restriction that all ants travel at a fixed speed, call it $v$. This assumption is reasonable because there is little variance in the cruising speeds of ants. Adopting it also allows us to define $\rho$ as a function simply of the angle $\theta$ of its velocity rather than of the vector $\vec{v}$, so that $\rho=\rho(\vec{x},\theta,t)$.

To define an equation for the evolution of $\rho(\vec{x},\theta,t)$, we must consider the factors that might change the density of particles at a particular location $\vec{x}$ and moving in the specific direction $\theta$ after an infinitesimal step in time. We draw this from the taxis model described by Codling, et. al. in their 2008 review paper on random walk models in biology. In their model, there exists some probability per time $\alpha$ that the velocity $\vec{v}$ of a particle will change at any time. Further, a reorientation kernel $T(\vec{v},\vec{v}^{\prime})$ is defined to be the probability that, given a reorientation occurs, the velocity of the particle will change from $\vec{v}^{\prime}$ to $\vec{v}$. Since $T(\vec{v},\vec{v}^{\prime})$ gives a probability, it satisfies the normalization condition that $\int T(\vec{v},\vec{v}^{\prime})=1$. The evolution equation arising from these definitions in the Codling paper is

[TABLE]

where $\nabla_{x}$ denotes the spatial gradient operator. Converting to our form of the system with fixed speed $v$ and $\rho(\vec{x},\theta,t)$ instead of $\rho(\vec{x},\vec{v},t)$, we have

[TABLE]

In this equation, the left-hand side is the total derivative of $\rho$ with respect to time, equal to the partial derivative of $\rho$ in time plus the convective derivative ($\vec{v}\cdot\nabla_{x}\rho$) in the center expression. The right-hand side of the equation defines the mechanism of change in the particle density $\rho$: The first term on the right side, $-\alpha\rho$, signifies that $\alpha\rho$ is the rate at which particles at position $\vec{x}$ and traveling in direction $\theta$ randomly change direction, and hence are no longer part of the $\rho(\vec{x},\theta,t)$ phase density. The second integral term can be understood as the rate at which particles traveling in directions $\theta^{\prime}\neq\theta$ randomly change their direction to $\theta$. Thus, we have a simple form for an evolution equation for the particle density $\rho$ in this definition.

The integration of both memory and reinforcement into our system, which we believe will produce the spiraling solution of the ant mill, is made in the specific form of the reorientation kernel $T$. In the Codling paper, the kernel $T$ is defined only in terms of the original direction $\theta^{\prime}$ and the new direction $\theta$. Clearly, this setup serves as a strong model for a correlated random walk, as the probability of a particle engaging in that direction change could be made higher for small values of $|\theta-\theta^{\prime}|$ and lower for large values of $|\theta-\theta^{\prime}|$. We seek, however, a kernel $T$ that combines correlation with reinforcement based on the surrounding concentrations of chemical. In order to construct such a kernel, we must first define an evolution equation for the chemical concentration. We give a simple chemical evolution equation below, where $g(\vec{x},t)$ denotes the concentration:

[TABLE]

In this equation, in which $\gamma$ and $\beta$ are constants, the first integral term indicates the rate at which the particles located at $\vec{x}$ deposit chemical, and the second term signifies the rate at which the chemical evaporates. In essence, the equation describes that the rate at which the chemical concentration is changing at a location $\vec{x}$ at a time $t$ is equal to rate at which ants are adding chemical minus that at which the chemical is evaporating (an exponential rate). Given this definition, we introduce the following form for the reorientation kernel $T$:

[TABLE]

where $\theta_{g}$ is the angle of the $\nabla g$ vector and $J$ is a constant. The reasoning behind this choice of kernel is as follows: The closer to $\theta-\theta_{g}$ is to [math], the greater the value of $T(\theta,\theta_{g})$ due to the nature of the cosine function. This means that the closer the angle $\theta$, the intended level of direction, is to the gradient, the more likely it is for a particle initially in the $\theta^{\prime}$ direction to change direction to $\theta$. Although we could also used $\theta^{\prime}$ in our definition of $T$, using only $\theta$ and $\theta_{g}$ simplifies the kernel while still making sense physically, so we proceed with the current form. The constant $\frac{1}{4\pi^{2}}$ serves as a normalization factor, and $|J|<1$ so that $T>0$.

We hence have the following integral equation for $\rho$:

[TABLE]

Our next step is to attempt to solve this equation for a steady-state, radially symmetric solution that resembles the ant mill. The symmetry of this solution form makes it valuable to convert it from rectangular to polar coordinates $(r,\phi)$. We notice that in the steady-state solution, the direction of the chemical gradient $\theta_{g}$ is always equal to the polar angle $\phi$, because the higher concentration of chemical in the center of the spiral pulls the ants around in a circle. Using these facts, letting $\phi=\theta_{g}=0$ without loss of generality, and simplifying, we arrive at the following equation from our original model,

[TABLE]

where $P(ik,\theta)$ is the Fourier transform of $\rho$ (i.e. $P(ik,\theta)=\mathcal{F}(\rho)$).

The above integral equation is a homogeneous Fredholm equation of the second kind. According to Arkfen in Section $16.3$ of his book Mathematical Methods for Physicists, however, this equation has no nontrivial solution: For all space, [math] is a solution to the equation, so we cannot find a nonzero solution. This result means that the current formulation of our system of equations modeling ant motion has no steady-state. In hindsight, this was to have been expected—it was unlikely that a spiraling time-independent solution could have been obtained without a nonlinearity in the equations causing the ants to remain in orbit.

In addition, the inability of the integral equation model discussed in this section to produce a nontrivial steady-state solution suggests that the fundamental assumption underpinning the model—that the density $\rho$ is a function of both position $\vec{x}$ and velocity $\vec{v}$—is ill-fated. Indeed, storing a density value for every position and velocity in space is a lot of information, as we are basically keeping track of all of the individual particles that are moving around together in the domain. In the next section, we remedy this potential problem of too much information by developing a slightly simpler model. In our new model, we consider the average velocity of particles at every point in space, rather than all the velocities of the different particles at that point. Furthermore, the inability of our linear integral equation to produce a nontrivial solution tells us that nonlinearity plays a critical role in describing circular motion. Thus, this aspect of our model must also be fixed, by the addition of nonlinearity. Since nonlinear integral equations are difficult to work with, we turn to non-linear diffusion-advection equations to model our system instead, because they preserve much of the core physics while also being more pliable. In combination, these changes ultimately end up producing the steady-state result that we seek.

Based on what we have learned from the integral equation model, we now develop a new, nonlinear model based on diffusion-advection equations.

3 Diffusion-advection model

3.1 Construction

To develop our new model of ant motion, we draw from past models in which the density $\rho$ is a function of only position and time, not velocity. The basis of our new system is the following pair of pure reinforcement equations for the particle density $\rho$ and chemical concentration $g$ in a one-dimensional environment from Othmer and Stevens (1997):

[TABLE]

[TABLE]

The first equation describes the evolution of the particle density and the second describes the evolution of the chemical concentration. The constant $\lambda$ indicates the magnitude of the rate of chemical deposition by the ants, and $D$ is the diffusion constant. We chose this particular system to serve as the basis for our new model because Othmer and Stevens have shown analytically and numerically that simple solutions of these equations are asymptotically stable. Further information on these simple solutions and their stability results can be found in Other and Stevens’ 1997 paper. For now, we press onward to transform this pure-reinforcement system into a mixed memory-reinforcement one by adding a formulation of velocity.

In Othmer and Stevens’ model, there is no concept of velocity, as indicated by the absence of a velocity variable $\vec{v}$ in the equations. Essentially, particles immediately forget their velocity at every instant, so there is no memory in the system. We change this by first adding a convective derivative to the left-hand side of the $\rho$ equation, so the left-hand side becomes the total derivative rather than the partial derivative. We also, of course, convert the one-dimensional equation into a two-dimensional one by changing the $\frac{\partial}{\partial x}$ operators to $\nabla$ operators:

[TABLE]

In addition, we must define a velocity evolution equation to complete the memory-reinforcement system. It seems reasonable that the rate at which the velocity at a location $\vec{x}$ changes should be related to the concentration of chemical at that location. In particular, it would make sense that the derivative of the velocity $\vec{v}$ is in the same direction as the chemical gradient $\vec{g}$, since particles are attracted to higher concentrations of chemical, and that the steeper the gradient, the higher the rate of change. Taking these assumptions into account, we develop the simple velocity equation

[TABLE]

where $b$ is a constant.

Hence, we have constructed the following system of three partial differential equations that models memory-reinforcement ant motion from the original system by Othmer and Stevens:

[TABLE]

[TABLE]

[TABLE]

Our next step is to solve our new system for a steady-state solution that resembles the ant mill.

3.2 Steady-state solution

We seek a solution to our system that resembles the ant death spiral. Such a solution must be radially symmetric; in other words, the density $\rho$, the chemical concentration $g$, and the velocity $\vec{v}$ must have no $\theta$ dependence. Since each of $\rho$, $g$, and $\vec{v}$ must therefore only depend on the radial distance $r$, we can rewrite the system of equations 11, 12, and 13 as

[TABLE]

[TABLE]

[TABLE]

using the polar forms of the gradient, divergence, and material derivative operators, and where $v_{\theta}$ denotes the $\theta$ component of $\vec{v}$. Further noting that $\frac{d\rho}{dt}=0$, $\frac{\partial g}{\partial t}=0$, and $\frac{\partial\vec{v}}{\partial t}=0$ in the time-independent case, we rewrite the equations as

[TABLE]

[TABLE]

[TABLE]

We can rearrange equation 18 to get

[TABLE]

and we substitute this expression for $g$ into equation 17 to obtain the single-unknown equation

[TABLE]

To solve equation 21 for $\rho$, we first integrate both sides with respect to $r$, so we obtain

[TABLE]

[TABLE]

where $C_{1}$ is a constant. We multiply both sides of the equation by $(\alpha+\beta\lambda\rho)$ and simplify as follows:

[TABLE]

[TABLE]

[TABLE]

where the last equation follows by changing the sign of the arbitrary constant $C_{1}$. This final form of the equation is straightforward to solve for $\rho$ using separation of variables. We subtract $C_{1}\frac{\beta}{\alpha}\lambda\rho$ from both sides of the equation, and then divide the right-hand side by $(-C_{1}-C_{1}\frac{\beta}{\alpha}\lambda\rho)$ and $\frac{r}{dr}$ so that each side of the equation is in terms of one variable:

[TABLE]

[TABLE]

We then integrate both sides and simplify to obtain

[TABLE]

[TABLE]

[TABLE]

where $C_{2}$ is a constant. Equation 25 can be rearranged to obtain the final form

[TABLE]

where the constant $C_{2}$ is redefined as $e^{C_{1}C_{2}}$, where the $C_{2}$ in the exponent is the original $C_{2}$.

This solution for $\rho$ is reasonable because, although $\rho$ approaches infinity at $r=0$, the integral of $\rho$ converges for $r=1$ to $r=\infty$ for values of the exponent. This models the death spiral because it would not make sense for moving particles to be circling extremely fast in the center of a physical spiral. This is observed in real ant death spirals, in which there is a certain radius in the center of the spiral within which no ants appear to be circling (see the ant death spiral with the empty circular center in Figure $1$ and the schematic in Figure $2$).

From our solution for $\rho$, since $g=\lambda\rho$, we also have

[TABLE]

as the solution for $g$. We then use this $g$ solution to obtain the solution for $\vec{v}$ by manipulating equation 19 with a little algebra, as below:

[TABLE]

[TABLE]

[TABLE]

Hence, the final steady-state, radially symmetric solution to the entire system is

[TABLE]

[TABLE]

[TABLE]

Since the solution approaches infinity as $r$ approaches [math], we note that the radial domain must be $r>r_{a}$ for a positive $r_{a}$. For simplicity, when performing stability analysis, we limit the domain further by specifying an outer boundary $r_{b}$, so we have a finite domain of $[r_{a},r_{b}]$.

By considering the signs of the constants in this steady-state solution, we know that the solution must decay to [math] as $r\rightarrow\infty$. We arrive at this conclusion by first noting that each of the constants $\alpha$, $\beta$, and $\lambda$ are positive, based on the nature of the phenomena we are modeling and, specifically, the fact that the ants are attracted to (not repelled by) the chemical. This means that the coefficients $\frac{\alpha}{\beta\lambda}$ and $\frac{\alpha}{\beta}$ in the $\rho$ and $g$ equations are positive, so the coefficient $C_{2}$ of $r$ must also be positive for the density and chemical concentration to be positive. Since $C_{2}$ is positive, $C_{1}+\frac{\alpha}{\beta\lambda}$ must also be positive, because otherwise the velocities $\vec{v}$ will be imaginary due to the square root in the coefficient $\sqrt{bC_{2}(C_{1}+\frac{\alpha}{\beta\lambda})}$ in the $\vec{v}$ equation. Finally, since the $C_{1}+\frac{\alpha}{\beta\lambda}$ is positive, the powers of $r$ in all three equations must be negative, indicating that the solutions all decay as $r\rightarrow\infty$. The fact that these coefficients confirm our expectation of what the solution should be is a promising sign of the robustness of the model.

Having arrived at the steady state, we can now perturb these time-independent solutions slightly and observe how the system behaves over time in response to determine its stability.

3.3 Stability Analysis

Before beginning stability analysis, we shall briefly review what stability analysis is and what it can show about a system. Essentially, in a time-dependent system of differential equations, a solution to the system can either blow up, collapse, or oscillate as time approaches infinity. Blowup refers to the scenario when the solution approaches positive or negative infinity as time increases. Collapse indicates when a solution approaches a finite value as time increases. In a system that oscillates as time approaches infinity, the limit of the solution does not exist (for $t\rightarrow\infty$).

To analyze the stability of our system of partial differential equations around the obtained steady-state solutions, we first let the steady-state solutions for $\rho$, $g$, and $\vec{v}$ be represented by $\rho_{0}(r)$, $g_{0}(r)$, and $\vec{v}_{0}(r)$. We define $\tilde{p}(r,\theta,t)$, $\tilde{g}(r,\theta,t)$, and $\tilde{v}(r,\theta,t)$ to be perturbations to the steady state equations, each of which is a function of $r$, $\theta$, and $t$. Each perturbation is of order $\epsilon_{1}>0$ where $\epsilon_{1}$ is a small real number. Thus we can write

[TABLE]

[TABLE]

[TABLE]

We seek expressions for the unknown perturbations $\tilde{p}$, $\tilde{g}$, and $\tilde{v}$ that solve the general (not time-independent and axi-symmetric) system. Analyzing these expressions will allow us to assess whether small perturbations to the steady-state solutions cause the system to blow up or stabilize in infinite time. To solve for $\tilde{p}$, $\tilde{g}$, and $\tilde{v}$, we will first linearize the time-dependent equations for $\rho$ and $\vec{v}$, both of which are nonlinear. We do this because it is far more straightforward to solve linear differential equations than non-linear differential equations, both analytically and numerically. The full details of the linearization process for the $\rho$ and $\vec{v}$ equations are given in Appendix A; in this section, we omit some of the algebra for the sake of concision.

First, we convert the time-dependent equation for $\rho$ to polar coordinates, a more convenient coordinate system to use for modeling the spiraling motion of ants:

[TABLE]

[TABLE]

Note that part of the nonlinearity in the equation 34 comes from the expression $(\frac{\beta}{\alpha+\beta g})$. We eliminate this nonlinearity by linearizing $g$ around the solution $g_{0}$, as follows:

[TABLE]

[TABLE]

Since $\epsilon_{1}$ is small, we can ignore terms of order $\epsilon_{1}^{2}$ and greater, so we have

[TABLE]

Substituting equation 36 into equation 34, along with the three perturbation equations $\rho=\rho_{0}+\epsilon_{1}\tilde{\rho}$, $g=g_{0}+\epsilon_{1}\tilde{g}$, and $v=v_{0}+\epsilon_{1}\tilde{v}$, we obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We notice that all terms in the above equation without an $\epsilon_{1}$ coefficient drop out, because $\rho_{0}$, $g_{0}$, and $\vec{v}_{0}$ form a time-independent solution ($\frac{\partial\rho_{0}}{\partial t}=0$). In addition, all derivatives of $\rho_{0}$, $g_{0}$, and $\vec{v}_{0}$ in $\theta$ are zero, because the radially symmetric solution has no $\theta$ dependence. Finally, all terms with a coefficient of $\epsilon_{1}^{2}$ can also be dropped since $\epsilon_{1}$ is small and thus terms that are second order in $\epsilon_{1}$ are negligible. Simplifying as such, and then dividing out the $\epsilon_{1}$ terms and rearranging, we get

[TABLE]

[TABLE]

[TABLE]

We can write a similar linear partial differential equation for $\frac{\partial g}{\partial t}$ by plugging $g=g_{0}+\epsilon_{1}\tilde{g}$ and $\rho=\rho_{0}+\epsilon_{1}\tilde{\rho}$ into the $\frac{\partial g}{\partial t}$ equation and canceling the non-$\epsilon_{1}$ coefficient terms, as below:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We linearize the $\vec{v}$ equation by a similar method, given explicitly in Appendix A.$2$, and we obtain the following two equations for $\tilde{v}_{r}$ and $\tilde{v}_{\theta}$:

[TABLE]

[TABLE]

from the $\hat{r}$ and $\hat{\theta}$ components of $\vec{v}$.

Thus, our linearized linear system of partial differential equations is the following:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

To prove asymptotic stability, we assume solutions of the form

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $s$ and $n$ are arbitrary parameters and $F$, $G$, $H_{r}$, and $H_{\theta}$ are functions of $r$. It is reasonable to assume that the perturbations in $\tilde{\rho}$, $\tilde{g}$, and $\tilde{v}$ have the same time and $\theta$ dependence—this assumption also greatly simplifies the upcoming analysis. We choose these particular solution forms because they give rise to the following useful properties:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

When substituted into the linear system of partial differential equations in terms of $\tilde{\rho}$, $\tilde{g}$, and $\tilde{v}$ and their derivatives in $r$, $\theta$, and $t$ (equations 42 to 45), these properties produce a system of ordinary differential equations with derivatives only in $r$. This system of ordinary differential equations can then be solved numerically for certain values of the parameters $s$ and $n$. If solutions satisfying the boundary conditions exist only for $s$ such that $s<0$, then we have shown that the system is asymptotically stable in time. If solutions exist only for $s$ such that $s>0$, then the system blows up in infinite time.

We substitute the properties into the system of partial differential equations to obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Equation 49 can be rearranged to get

[TABLE]

Substituting this expression for $\tilde{g}$ into equation 48, we simplify so that the equation is only in terms of $\tilde{\rho}$, $\tilde{v}$, and their derivatives:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Finally, with a little algebra, we can manipulate equations 49, 50, and 51 to get the following expressions for $\tilde{v}_{\theta}$ and $\tilde{v}_{r}$:

[TABLE]

[TABLE]

Note that the expression for $\tilde{v}_{\theta}$ is only in terms of $\tilde{v}_{r}$, and the expression for $\tilde{v}_{r}$ is only in terms of $\tilde{\rho}$. Since our linear equation for $\tilde{\rho}$ gives an expression in terms of only $\tilde{\rho}$ and $\tilde{v}_{r}$, substituting equation 55 in the $\tilde{\rho}$ equation gives us a single ordinary differential equation in terms of only $\tilde{\rho}$ and its derivatives.

We have hence obtained a system of ordinary differential equations from our original system of partial differential equations. This system can be solved numerically, and the obtained numerical solution can then be analyzed to prove stability or instability. We explain our numerical methods below.

We divide the domain $[r_{a},r_{b}]$ of our system into $N-1$ spaces using $N$ division points (including the endpoints $r_{a}$ and $r_{b}$). We label these division points $r_{1}$, $r_{2}$, $\ldots$, $r_{N}$. We define a time step $\Delta t$, so that the total elapsed time $t=\Delta t\cdot i$ for each index $i$. For each time $t=\Delta t\cdot i$, we define a vector $\vec{w}_{i}$ of size $4N\times 1$ as follows:

[TABLE]

Note that $r_{1}=r_{a}$ and $r_{N}=r_{b}$ in the above vector.

We now discretize in space, so we get

[TABLE]

where $M$ is a $4N\times 4N$ block matrix with $16$ $N\times N$ blocks:

[TABLE]

Each of the $N\times N$ blocks corresponds to a coefficient of one of the terms in the linear system of equations (i.e. either a coefficient of $\tilde{\rho}$, $\tilde{g}$, $\tilde{v}_{r}$, and $\tilde{v}_{\theta}$).

Once all the matrices are filled, we run code that constructs the matrix $I-M\Delta t$ for a fixed parameter and spacing $\Delta r$, and then calculates the norm of the matrix $I-M\Delta t$. If the norm is greater than $1$, then the system is stable. We found that for large values of the constant $b$, this was indeed the case, indicating that there exist conditions under which our steady-state solution is stable. Thus, we have shown numerically that the naturally-arising ant mill is a stable configuration.

4 Conclusions

In this paper, we have developed a system of diffusion-advection partial differential equations that models the spiraling configuration of ants known as the ant mill. Our system integrates both the biases of memory (correlation), in which moving particles remember their past velocities and tend to maintain them, and positive reinforcement, in which particles travel up chemical gradients. We have shown analytically that there exists a steady-state solution to our system that resembles the ant mill, and have proven numerically that this steady-state solution is stable. Moreover, we have also developed and discussed the inadequacies of integral model equations to describe mixed memory- and reinforcement-based motion.

The creation of a continuous random walk model that mixes memory and reinforcement, shows potential to model real-world memory-reinforcement systems, and which lends itself nicely to stability analysis is an important result in random walk theory. Previously, most random walk models of natural phenomena were based either purely on memory or purely on reinforcement. This limitation caused the mathematical models developed to describe many phenomena to be slightly or significantly mismatched with the observed systems. Now, given a simple mixed memory-reinforcement model, further extensions can be made that allow a much larger number of random walk scenarios to be mathematically described.

The applications of a continuous mixed memory-reinforcement model go far beyond modeling the spiraling motions of ants. Several species of bacteria, such as the soil bacteria known as myxobacteria, as well as slime molds also move by a combination of memory-based and reinforcement-based motion. Furthermore, and perhaps most excitingly, memory and reinforcement also jointly influence the growth and directions of movement of new blood vessels in the process of angiogenesis. This process is significant because it occurs during the growth of a tumor in an organism, as the cancerous cells of the tumor attract growing blood vessels to extend toward them to provide nourishment. If the specific values of the parameters in the memory-reinforcement model that replicates observed patterns of angiogenesis are determined, we will have understood the conditions of angiogenesis mathematically. Such an understanding could lead to methods of manipulating the growth of blood vessels around tumors so that they do not grow toward the cancer cells, potentially resulting in new cancer-fighting treatments.

In addition to its potential in developing cancer treatments, the mixed memory and reinforcement model developed in this paper also can be applied to the more humble purpose of determining the conditions under which ants can become locked into a death spiral, and the probability of such a configuration occurring. The specifications and extensions of our prototype model to describe all of these different natural systems are what makes it a relevant and interesting topic of study, and are left to future work.

Appendix A Linearization of Time-Dependent $\rho$** and $\vec{v}$** Equations

A.1 Density Equation, $\rho$

We begin with the density equation $\rho$ in Section 3 based on the 1997 paper by Othmer and Stevens,

[TABLE]

where $\rho$ represents the particle density, $g$ represents the pheromone concentration, and $\vec{v}$ the velocity of the particles. We expand the divergence operator on the right-hand side of the equation to obtain

[TABLE]

Next, we convert the equation to cylindrical coordinates with radial variable $r$ and angle variable $\theta$:

[TABLE]

[TABLE]

Taking the derivatives on the right-hand side and simplifying, we get the following equation:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We now linearize the obtained nonlinear equation. We introduce the notation $\rho_{0}$, $g_{0}$, and $\vec{v}_{0}$ to denote the steady-state solutions of the system of equations defined in Section 3, and $\tilde{\rho}$, $\tilde{g}$, and $\tilde{v}$ to denote perturbations to the steady-state. The perturbations are of order $\epsilon_{1}>0$ for small $\epsilon_{1}$, so we have

[TABLE]

[TABLE]

[TABLE]

We also note the nonlinearity in the above equation due to the expression $\frac{\beta}{\alpha+\beta g}$. We resolve this issue by linearizing $g$ around the solution $g_{0}$ so that

[TABLE]

[TABLE]

Since $\epsilon_{1}$ is small, we can ignore terms of order $\epsilon_{1}^{2}$ and greater, so we have

[TABLE]

Substituting this and equations 63 into the nonlinear equation, we obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We now note that all terms without an $\epsilon_{1}$ coefficient drop out of the equation because $\rho_{0}$, $g_{0}$, and $\vec{v}_{0}$ form a time-independent solution ($\frac{\partial\rho_{0}}{\partial t}=0$). We also note that all derivatives of $\rho_{0}$, $g_{0}$, and $\vec{v}_{0}$ in $\theta$ are zero because the axi-symmetric solution has no $\theta$ dependence. Finally, all terms with a coefficient of $\epsilon_{1}^{2}$ can also be dropped since $\epsilon_{1}$ is small. Simplifying as such, we get

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Dividing out the $\epsilon_{1}$ terms, we therefore have

[TABLE]

[TABLE]

[TABLE]

which is linear. We can rearrange this equation so that all terms with $\tilde{\rho}$, $\tilde{g}$, $\tilde{v}$, and their derivatives on the right-hand side are grouped to obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We now assume solutions of the form

[TABLE]

[TABLE]

[TABLE]

[TABLE]

which we will use to prove asymptotic stability. The variables $s$ and $n$ are arbitrary parameters and $F$, $G$, $H_{r}$, and $H_{\theta}$ are functions of $r$. These solutions give rise to the following useful properties:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Plugging in these properties into our equation, we get

[TABLE]

[TABLE]

[TABLE]

[TABLE]

which we can rearrange with a little algebra to obtain

[TABLE]

[TABLE]

[TABLE]

Finally, we plug the relation $\tilde{g}=\left(\frac{\lambda}{s+1}\right)\tilde{\rho}$ from Section 3 (equation 52) and rearrange to get an equation only in terms of $\tilde{\rho}$ and its $r$ derivatives:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A.2 Velocity Equation, $\vec{v}$

We start with the velocity equation $\vec{v}$ defined in Section 3:

[TABLE]

Taking the material derivative $(\vec{v}\cdot\nabla)\vec{v}$, we obtain

[TABLE]

which we rearrange so that all $\hat{r}$ and $\hat{\theta}$ terms are grouped to get

[TABLE]

[TABLE]

Separating the $\hat{r}$ and $\hat{\theta}$ components, we have the following two equations:

[TABLE]

[TABLE]

We begin by linearizing equation 78.

Substituting $v_{r}=v_{0,r}+\epsilon_{1}\tilde{v}_{r}$ and $v_{\theta}=v_{0,\theta}+\epsilon_{1}\tilde{v}_{\theta}$, we obtain

[TABLE]

[TABLE]

Noting that $v_{0,r}=0$ and eliminating all terms that are zeroth or second order in $\epsilon_{1}$ (since $\frac{\partial\vec{v}_{0}}{\partial t}=0$ and $\epsilon_{1}$ is small), we obtain

[TABLE]

Dividing out $\epsilon_{1}$, we get

[TABLE]

We now substitute the solutions for $\tilde{v}_{r}$ and $\tilde{v}_{\theta}$ in $(1)$ into the equation to obtain

[TABLE]

which we rearrange to get

[TABLE]

We now perform the same manipulations on equation 79. Substituting $v_{r}=v_{0,r}+\epsilon_{1}\tilde{v}_{r}$ and $v_{\theta}=v_{0,\theta}+\epsilon_{1}\tilde{v}_{\theta}$, we obtain

[TABLE]

[TABLE]

Since $v_{0,r}=0$ and we eliminate all terms that are zeroth or second order in $\epsilon_{1}$, we get

[TABLE]

Dividing out $\epsilon_{1}$ gives us

[TABLE]

Finally, plugging in the solutions for $\tilde{v}_{r}$ and $\tilde{v}_{\theta}$ as above and rearranging gives us

[TABLE]

[TABLE]

We therefore have the following system of linear equations in $\tilde{v}_{r}$, $\tilde{v}_{\theta}$, and $\tilde{g}$:

[TABLE]

[TABLE]

We solve this system and plug in the relation $\tilde{g}=\left(\frac{\lambda}{s+1}\right)\tilde{\rho}$ (equation 52) to get

[TABLE]

[TABLE]

which we used in Section $3.3$.