Analytical solution of the systems of nonlinear fractional partial differential equations using conformable Laplace transform iterative method

1.   Introduction and background

P. polycephalum, or a slime mold, is an amoeboid organism that is notable for its ability to perform complex tasks despite its relatively simple biological structure and lack of a brain. Surprisingly, slime molds are known to be able to solve mazes, construct robust networks, and solve shortest path and spanning tree problems ^{[1,2,3,4,5,6,7,8,9]}. The simplicity of the slime mold organism and the complexity of its emergent behaviors have inspired research into P. polycephalum as a model for collective behavior ^{[10,11,12,13]}. A hallmark of research in collective dynamics is complex global structure emerging as a result of relatively simple local interactions between agents. Indeed, cellular automata models inspired by P. polycephalum food seeking behavior have been found to be good models of other collective phenomena such as the formation of transportation networks in a country or the "cosmic web" of stellar material between galaxies ^[14,15]. When seeking food, P. polycephalum moves to concentrate most of its mass on food sources - however it maintains a network of mass between food sources. In the absence of food it begins to retract to a central mass. Past research into models of P. polycephalum has mainly attempted to model the general food seeking behavior of physarum or its food seeking behavior within a maze. In the first camp most of the models take a microscopic, cellular automata approach - the slime mold is represented by a collection of particles that evolve discretely in time, space, or both. Different interaction rules and sometimes even different "types" of agents are employed to model different slime mold behaviors, however almost all models include a rule that causes slime mold agents to move towards food sources ^{[16,17,18,19,20]}. Most attempts to model the maze solving behavior of the slime mold take a macroscopic approach. Some studies note that the P. polycephalum moves and forms networks by moving fluid through tubules within its mass - much like an electrical network ^[10,21,22]. These authors focus on modeling P. polycephalum as a network of ordinary differential equations that describe the flow of P. polycephalum through network components. Others try to be more specific and attempt to directly model the density of P. polycephalum on each maze component. Here the maze is described as a network topology and P. polycephalum is characterized with one dimensional partial differential equations for chemotaxis (in particular the Keller-Segel model) defined on the edges of the network with coupling conditions at the junctions of the maze ^[23,24]. Both of these macroscopic descriptions attempt to model the path that P. polycephalum builds between food sources once they are located. In particular they do not attempt to model the general food seeking behavior of the organism.

In this manuscript we propose a model hierarchy aimed at modeling the food seeking behavior of P. polycephalum. The main behaviors we attempt to model are agammaegation of mass on food sources while maintaining connected paths between food sources and retraction in the absence of food sources. We first define an agent based model that is continuous in time and space through a system of stochastic differential equations. The evolution of each agent is governed by three terms-an interaction term with other agents that causes agents to agammaegate, a drift term that causes agents to move towards food sources and a noise term intended to model slime mold foraging behavior in the absence of food sources. However any simulation of the microscopic model that includes a suitable number of agents to accurately portray P. polycephalum behavior would be computationally intensive. Therefore, we demonstrate under the assumption of propagation of chaos that as the number of agents approaches infinity, the marginal distribution of each agent converges to a deterministic density based model in the form of an agammaegation-diffusion equation. This macroscopic model includes terms that correspond to the drift, interaction and noise terms present in the microscopic model. We then embark on an analysis of the macroscopic equation from the point of view of modeling P. polycephalum food seeking behavior. We first discuss a framework for choosing an appropriate interaction kernel through an analysis of the steady states of the macroscopic model in the case that there are no food sources. Then, we rigorously show that there is a parameter regime in which the only possible steady state of the macroscopic model is zero and thus this regime is not suitable for modeling P. polycephalum food seeking behavior. Finally, we simulate the macroscopic model using the blob method for agammaegation-diffusion equations presented in ^[25] and find that it is necessary to scale the three main terms of the model in order to reproduce agammaegation on food sources. We then simulate the evolution of the model in different scaling regimes and with different choices of interaction kernel and qualitatively examine which regimes result in the main behaviors we attempt to model. We find that there are several regimes that result in the qualitative desired behaviors.

From a modeling point of view our approach distinguishes itself from the previously described literature in several ways. As previously discussed, past models of P. polycephalum behavior mainly took either a discrete agent based or macroscopic viewpoint - our approach unifies these strategies as our macroscopic description is formally derived from the agent based model. Additionally, most models of general P. polycephalum food seeking behavior were from the agent based view point; our model hierarchy allows us to present a macroscopic description of this behavior. Finally, the macroscopic description that results is an agammaegation diffusion equation with a drift term; while known, to our knowledge such models have not been widely studied. Agammaegation diffusion equations without a drift term are well known to exhibit "blow up" below a critical mass - we do not observe such a phenomena in our simulations (although we do not vary the mass) - it could be an interesting area of future research to investigate the effect of the drift term on the occurrence of blow up.

2.   Model derivation

The main features of slime mold food seeking behavior are agammaegation of mass on food sources while maintaining a connected mass and retraction in the absence of food sources. To model these behaviors we first introduce an agent-based model. Here, the slime mold is represented by a swarm of $N$ agents who's trajectories, $\mathbf{{X}}_{i}(t) \in \mathbb{R}^{d}$, are modeled as a system of stochastic differential equations. In the following we will refer to this model as the microscopic model in order to reinforce the intuition that the agents represent slime mold "particles".

Definition 1 (Model 1-Microscopic Model). Consider a collection of $N$ agents with (random) positions

$(\mathbf{{X}}_{1}(t), ..., \mathbf{{X}}_{n}(t))\in (\mathbb{R}^{d})^{N}$. The microscopic model is given by the evolution:

where $\{\mathbf{{B}}_{i}\}_{i}$ are independent, standard $\mathbb{R}^{d}$-valued Brownian motions, $V, W: \mathbb{R}^{d}\rightarrow \mathbb{R}$ are $C^1$ functions, $\{\mathbf{{X}}_{i}(0)\}_{i}$ are independent and identically distributed random variables independent of $\{\mathbf{{B}}_{i}\}$ with probability density function given by $\rho_{0}: \mathbb{R}^{d}\rightarrow \mathbb{R}$, and $A, B, C > 0$.

We note that if we assume that $\nabla V$ and $\nabla W$ are globally Lipschitz and further if there exists $K_{1}, K_{2} > 0$ such that for $x \in \mathbb{R}^{d}$:

then it can be shown that Eq (2.1) has a unique strong solution.

Intuitively, $V$ represents the density of chemo-attractants given off by food sources. For this reason, we usually think of $V$ as a sum of radially symmetric, positive functions. Each term represents the attractants emanating from a single food source. A prototypical example would be a sum of Gaussians centered at different positions. $W$ represents the interaction between agents. Notice that, roughly, if $\nabla W(\mathbf{{X}}_{j}(t) - \mathbf{{X}}_{i}(t)) \approx \mathbf{{X}}_{j}(t) - \mathbf{{X}}_{i}(t)$, then the interaction causes agents to attract each other. Likewise if $\nabla W(\mathbf{{X}}_{j}(t) - \mathbf{{X}}_{i}(t)) \approx \mathbf{{X}}_{i}(t) - \mathbf{{X}}_{j}(t)$, then the agents will repel each other. As we have noted, P. polycephalum tends to maintain a connected mass, even when agammaegating around disparate food sources. Therefore, this interaction term should cause agents to remain locally close, i.e. agents should not repel each other at long ranges. The noise term, $Cd\mathbf{{B}}_{i}(t)$, models P. polycephalum foraging behavior. In the absence of chemo-attractants, P. polycephalum forages by diffusing outward from its initial position while still maintaining a connected mass. Here, in the absence of a food source ($V = 0$) the noise term should cause particles to spread from their initial positions while the interaction term causes them to still have some propensity to remain locally "together".

We include the parameters $A, B$ and $C$ for the case that simulations of Eq (2.1) demonstrate the need to scale the terms' influence on each particles trajectory to produce evolutions that resemble P. polycephalum. However, any simulation of Eq (2.1) that would give an accurate picture of how Eq (2.1) models P. polycephalum would have to use a very large number of particles; this becomes computationally intractable. Therefore, we must find a way to examine the behavior or Eq (2.1), for a large number of particles, without computing the trajectory of every particle explicitly. We will show that in the large particle limit that the evolution of the marginal distribution of any particle is given by an agammaegation-diffusion equation. There exist computationally efficient methods for simulating such equations, therefore we can examine the behavior of Eq (2.1) by simulating the evolution of the marginal density. The technique we will employ to derive an appropriate model is the assumption of propagation of chaos; i.e. if the number of agents agents is large, the trajectories of any two agents can be assumed to be independent as they are not likely to interact strongly. This is true in many cases, see for example ^[26,27].

Theorem 1. In the $N\rightarrow \infty$ limit and for a given, $A, B, C > 0$, $W, V \in C^1(\mathbb{R}^{d})$ and $\rho_{0}$, the evolution of the marginal distribution of any agent evolving according to Eq (2.1) is a solution of the agammaegation-diffusion equation

under the assumption of propagation of chaos; that is for any finite group of $M$ agents there exist distributions $\rho_{1}, ..., \rho_{M}$ on $\mathbb{R}^{d}$ such that the joint distribution $\rho^{M}(x_{1}, ..., x_{M}, t)$ satisfies:

for any time.

Proof. Consider a finite group of $N$ agents, $(X_{i}(t), t\geq 0)\subseteq \mathbb{R}^{d}, i = 1, ..., N$ evolving according to Eq (2.1). Let $\mu^{N}$ be the empirical distribution of the configuration of agents. That is:

where $\delta_{X_{i}(t)}$ is the Dirac measure with unit mass at $X_{i}(t)$. Let $F$ be a test function in $C_{c}^{\infty}((\mathbb{R}^{d})^{N})$, therefore by definition of the Dirac measure we deduce:

we first derive an expression for the time evolution of the expectation of this quantity. Let $\mathbf{{X}}(t) = (X_{1}(t), ..., X_{N}(t))$. By conditioning we find

where $\rho^{N}(\mathbf{{x}}, t)$ is the joint distribution of the collection of agents. Using Ito's Lemma (e.g. ^[28]) we can further expand the conditioned term in the above

We first focus on the third term of Eq (2.3). Notice that as each $X_{i}$ evolves according to Eq (2.1), since the noise terms for disparate agents are independent, we have that:

Therefore, the third term of Eq (2.3) simplifies as:

Now, turning to the second term of Eq (2.3) we find by appyling Eq (2.1) that:

As the stochastic integral is a martingale, we have that the last term of Eq (2.5) is equal to $0$. As a result, by combining Eqs (2.4) and (2.5), we have that:

So, by Eq (2.2) and the continuity of the trajectories, Eq (2.6) implies $P$ almost-surely that:

Therefore, $\rho^{N}$ is a weak solution to the following initial value problem:

Using Eq (2.8), we now compute the time evolution of the marginal distribution of one agent in the $N\to \infty$ limit under the assumption of propagation of chaos, that is for any finite collection of $M$ agents we have that their joint distribution satisfies (in a limiting sense):

In the following we consider a finite group of $M$ particles in the $N\to\infty$ limit. Without loss of generality we will consider the marginal distribution, $\rho(x_{1}, t)$ of the first agent which is given by:

and therefore, by Eq (2.8), we have that $\rho(x_{1}, t)$ satisfies (weakly):

We analyze Eq (2.10) term by term. First, by applying Eq (2.9) we find that

Next, if we again apply Eq (2.9) and integrate by parts in each term of the sum where $j\neq1$ we find that

Similarly we deduce

We can simplify Eq (2.13) further by integrating under the divergence operator and applying Eq (2.9), notice that:

Hence:

Therefore, combining Eqs (2.11), (2.12) and (2.15) we have by Eq (2.10) that the marginal of any agent satisfies:

as desired.

As the equation derived in Theorem 1 describes how the entire collection of particles evolves from a density standpoint, in the following we will refer to it as the macroscopic model.

Definition 2 (Model 2 - Macroscopic Model). The macroscopic model is given by the evolution of the agammaegation-diffusion equation

for $W, V \in C^1(\mathbb{R}^{d})$ and $\rho_{0}\in H^1(\mathbb{R}^{d})$ with $\int_{\Omega} \rho_{0}(x)dx = 1.$

Here, $V$ and $W$ continue to represent the chemoattractant of food sources and the slime mold's propensity to agammaegate respectively. In the following we will study Eq (2.17) from the viewpoint of using it to model P. polycephalum food-seeking behavior.

3.   Model properties

Many questions remain about how to utilize the macroscopic equation derived in the previous section to model a slime mold. In this section we present some analysis of stationary states of the equation in order to inform some of these modeling choices. First, we investigate stationary states in the case where there is no food present as a way to gain some information about what choice of kernel function might be reasonable. Inspired by ^[2], we believe that Gaussians represent a reasonable model for stationary states of a slime mold - however we find that Gaussian stationary states are only possible in the case that the kernel function, $W$, is quadratic. Next, we investigate stationary states of the scaled equation. We find through a fixed point argument that if the diffusion scaling parameter, $C$, is sufficiently higher than $A$ and $B$ then that the only stationary state of the equation is $0$ (even in the presence of a food source). This suggests that when tuning the scaling parameters that there cannot be "too much" diffusion.

3.1.   Zero-food stationary states

In the case that there are no food sources a slime mold will "retract" to a more compact configuration ^[2,29,30]. Often, this configuration is roughly radially symmetric. Given that and inspired by the central limit theorem, we believe that a reasonable model for the configuration of the slime mold in the case that there is no food could be a Gaussian. We now investigate if there are conditions on the kernel function that are imposed by the assumption that stationary states are Gaussian. For simplicity we work in one dimension however the calculations are analogous in higher dimensions. In the case that there are no food sources Eq (2.17) becomes:

Therefore a stationary state satisfies (in the following we will ignore the scaling parameters, $A$ and $B$, as they do not change the computation):

which implies that

for some constant $D$. Finally we have that:

Note that Eq (3.2) represents a general condition that zero food stationary states must satisfy. If we assume that stationary states are Gaussian, i.e. that:

then by Eq (3.2) we have that:

The above can only hold if $W'(y) = ay$ for some $a\in \mathbb{R}$. Therefore, Gaussians are zero food stationary states only if the interaction kernel is quadratic. Further, by plugging in $ay$ for $W'(y)$ in the above we can see that the following must hold:

So, for a given quadratic kernel there are a family of Gaussians that are potential zero food stationary states (depending on the mass of the initial profile). For an illustration of this fact see Figure 1. Here, we simulate Eq (2.17) in the case of no drift term using the "blob method" for agammaegation-diffusion equations introduced in ^[25]. We will discuss the numerical method in more detail in the Numerics section of the paper. We choose an initial condition representative of a slime mold agammaegated around two food sources at $x = 1$ and $x = -1$. We find in accordance with our calculation above that the mass profile appears to converge to a Gaussian configuration. This also illustrates a modeling property of Eq (2.17); in the case that food sources "run out" the profile will "retract" to a central configuration.

Effectively, by starting from information about the zero food stationary state, we have "solved" for the correct kernel. This analysis suggests an empirical method for determining the correct kernel function from experimental data. Instead of our assumption that zero food stationary states are Gaussian, statistical analysis of slime mold configurations in the presence of zero food sources could provide information about the stationary state configuration. This information could then be used to "solve" for the kernel that results in this "correct" stationary state.

3.2.   High diffusion stationary states

We now turn to considering stationary states of the scaled equation:

We will see in a later section that the "fair" regime $(A = B = C)$ does not result in food-source agammaegation. Therefore we are motivated to study different scalings of the terms in Eq (2.17). Here, we will rule out a large regime of scaling parameters by employing a fixed point argument to show that if the diffusion parameter, $C$, is sufficiently larger then the drift and interaction parameters then the only stationary state of Eq (2.17) is $0$ on bounded domains.

Theorem 2. Let $\Omega\subseteq \mathbb{R}^{d}$ be compact and connected. For a given $A, B > 0$, $V, W \in H^{2}(\Omega)\cap L^{\infty}(\Omega)$ and $\rho_{0}\in H_{0}^{1}(\Omega)$ then for any sufficiently large $C > 0$, $0$ is the unique stationary state of:

in $H_{0}^{1}(\Omega)$.

Proof. Recall that $\rho$ is a stationary state of Eq (2.17) if:

We define:

Given $\rho \in H_{0}^{1}(\Omega)$ we know that $G(\rho) \in L^{2}(\Omega)$ and that the equation:

has a unique solution in $H_{0}^{1}(\Omega)$. Given $K > 0$, define $H_{0}^{1}(\Omega, K) : = \{\rho \in H_{0}^{1}(\Omega) : ||\rho||_{H_{0}^{1}} \leq K\}$. It can be shown (see, for example, Theorem 4 in section 6.3 of ^[31]) that if $\tilde{\rho}$ is the unique solution of Eq (3.4) then the following estimate holds for some $D > 0$:

Therefore, by applying the triangle inequality and Cauchy-Schwarz inequality we can see that for a sufficiently large choice of $C$ if $\rho \in H_{0}^{1}(\Omega, K)$ then we must also have that the solution to Eq $(3.4)$, $\tilde{\rho}$, belongs to $H_{0}^{1}(\Omega, K)$ as well.

Define $\tilde{A}:H_{0}^{1}(\Omega, K)\rightarrow H_{0}^{1}(\Omega, K)$ via:

We will show that $\tilde{A}$ is a contraction and thus has a unique fixed point by Banach's fixed point theorem. Since $0$ is trivially a fixed point of $\tilde{A}$ and $K$ is arbitrary this implies that $0$ is the only fixed point of $\tilde{A}$ and thus the only stationary state of Eq (2.17). Let $\tilde{\rho}_{1} = \tilde{A}[\rho_{1}]$ and $\tilde{\rho}_{2} = \tilde{A}[\rho_{2}]$. We will show that there exists a $0 < \gamma(A, B, C) < 1$ such that:

Since $\Omega$ is bounded by assumption we have by Poincare's inequality that:

for some $\alpha > 0$. Therefore since $\tilde{\rho_{1}}$ and $\tilde{\rho_{2}}$ are both solutions of Eq (3.4) we must have by another application of the Poincare inequality that:

Therefore we have that:

and by Eq (3.6) that:

We claim that

for some $\gamma > 0$ that depends only on $A, B$ and $C$ for a given $W$ and $V$. Expanding we find that:

We will continue our analysis term by term. Starting with the first term and applying the Cauchy-Schwarz inequality we find that

Now, turning to the second term of Eq (3.9) and applying the Cauchy-Schwarz inequality we find:

We will further expand the last term of Eq (3.11). Notice that:

Now, combining Eqs (3.11) and (3.12) we find

Finally, we turn to the third term of Eq (3.9). Notice that:

Now, combining our results in Eqs (3.10), (3.13) and (3.14) and applying the definition of the norm in $H_{0}^{1}(\Omega)$ we find that:

where $\gamma \geq 0$ and is independent of $\rho_{1}$ and $\rho_{2}$ as we have assumed a uniform bound on their norms (they are members of $H(\Omega, K)$). Therefore, if $W$ and $V$ are given, $\gamma$ can be made arbitrarily small if $C$ is sufficiently larger than $A$ and $B$. Additionally we have that

for some $\beta > 0$ and therefore

Combining Eqs (3.7) and (3.16) we obtain:

Therefore for $C$ sufficiently larger than $A$ and $B$, the mapping $\tilde{A}$ is a contraction as desired.

From a modeling point of view, Theorem 2 demonstrates that a sufficiently strong contribution from the diffusion term will destroy any of the qualitative behaviors we wish to capture. We will see this explicitly in our numerical studies of Eq (2.17) in the next section.

4.   Numerics - Qualitative features of the macroscopic model

We now present some numerical investigations of the macroscopic model derived in the previous section.

These investigations are aimed at replicating qualitative properties of slime mold movement seen in ^[2,29,30], specifically agammaegation around food sources while still maintaining a "connected" mass. Recall that here, the function $V$ models the density of chemoattractants dispersed by food sources. Therefore, $V$ will always take the form of a sum of radially symmetric positive functions where each term is centered on a food source. We simulate Eq (4.1) via the "blob method for agammaegation-diffusion equations" introduced in ^[25] which approximates Eq (4.1) by solving the ODE:

for a collection of $N$ particles, $(x_{1}(t), ..., x_{N}(t), t\geq 0)$, who's initial positions are a regular grid on the domain on which we'd like to approximate Eq (4.1) and who's masses are given by $m_{i} = \rho(0, X_{i}(0))$. To recover an approximation of $\rho$ from the positions of the particles we convolve the particle solution with a mollifier ${\rm{ \mathsf{ φ} }}_{\epsilon}$:

We assume mollifiers are always of the form:

where $d$ is the dimension of the domain. It is shown in ^[25] that under some regularity and growth conditions on $W$ and $V$ that $\bar{\rho}$ converges to the solution of Eq (4.1). As a first example we will consider Eq (4.1) in one dimension. We will take the food source to be

Therefore we can think of the food sources as being located at $x = -1$ and $x = 1$. For the parameters of the ODE we take:

● $N = 100$ particles initially equally space on the interval $[-2.1, 2.1]$, i.e. with separation $h = .042$

● $\epsilon = h^{.99}$ (Mollifier parameter)

● initial mass profile given by $\rho_{0}(x) = \frac{1}{\sqrt{2\pi\sigma^{2}}}e^-\frac{x^{2}}{2\sigma^{2}}$, $\sigma^{2} = 0.0625$ as in ^[25].

Note that as the variance of the initial profile is very small that, at the particle level, slime mold particles have a very low probability of being found outside of $[-2.1, 2.1]$ (initially). In Figure 2 we perform a first simulation of Eq (2.17) in one dimension via the blob method in the so called "fair regime" ($A = B = C = 1$). Here, we choose a quadratic interaction kernel - specifically we impose that $W'(x) = x$.

Qualitatively, we do not observe agammaegation of slime mold mass around food sources to the degree present in ^[2]. Therefore, motivated by this example, in the following we will attempt to modify Eq (4.1) in order to better model slime mold food seeking behavior. In Figure 2 we observe that diffusion appears to be the dominating effect. Therefore, in order to examine the interplay of the diffusion term against the other two terms in Eq (4.1) we will employ the scaling parameters, $A, B, C > 0$, in order to control the contributions from the drift term, interaction term and diffusion term respectively. Roughly we examine two main regimes; "drift dominated" (where $A > B, C$) and "interaction dominated" (where $B > A, C$). We also examine so called "competition regimes" where two scaling parameters are equal but dominate the third. We do not examine the diffusion dominated regime as we've previously seen that if $A = B = C$ (we will refer to this as the "fair regime") then diffusion is the dominating effect. In each regime we examine three different interaction kernels - quadratic, polynomial and Gaussian. We choose the same ODE parameters as in the previous fair regime simulation.

4.1.   Food source dominated regime

We will first examine simulations of Eq (2.17) in the regime where the drift term (which models attraction to the food sources) is the dominating factor. We first examine the case where the kernel is quadratic and given by $W(x) = \frac{x^{2}}{2}$. Here, the kernel is purely attractive. At the agent based level all agents exert a "pull" on other agents with a force that is proportional to their spacial seperation. The results can be seen in Figure 3. Here, we qualitatively observe that in the case where $A = 10$, $B = C = 1$ that we have strong food source agammaegation however, the mass profile splits into two distinct bumps - the strong contribution from the drift term prevents a connected mass from being maintained. In the other three cases a connected mass is maintained and food source agammaegation is observed. However, in the two cases where $C = 5$ the degree of agammaegation is relatively weak likely due to the stronger contribution from the diffusion term.

Next, we examine the case where the kernel is polynomial and given by $W(x) = \frac{x^{4}}{4} - \frac{x^{2}}{2}$. Intuitively, at the particle level, this means that this kernel is attractive at long-ranges but includes a short-range repulsion to prevent particles from becoming too close. The results can be seen in Figure 4. Here, we observe results that are qualitatively very similarly to the results obtained in the quadratic case. This is not surprising given that the polynomial kernel is largely attractive and at long ranges the degree of attraction is proportional to the separation between particles - analogous to the quadratic kernel. However in the case where $A = 10, B = 5, C = 1$ the evolution with a polynomial kernel does not maintain a connected mass in contrast to the evolution witha quadratic kernel. This is likely due to the short-range repulsion that is present in the polynomial kernel. At the particle level, particles in the middle of the domain are repelled away from each otherwhile simultaneously being pulled by a food-source - encouraging a separation.

Finally, we examine the case where the kernel is Gaussian and given by $W(x) = {\rm{ \mathsf{ φ} }}_{.99}(x) = \frac{1}{(4\pi \cdot.99^{2})^{1/2}} e^{-x^{2}/4\cdot.99^{2}}$. Here, similarly to the quadratic kernel the interaction is purely attractive. However, unlike the polynomial the strength of attraction between two particles is not proportional to their separation. Instead, attraction becomes stronger (but not in an unbounded fashion) as two particles become closer and weaker as they become farther apart. This encodes the modeling assumption that interactions should be "local" - particles who are physically distant from each othershould not interact strongly. The results can be seen in Figure 5.

Here, we observe food source agammaegation in every case. similarly to the polynomial and quadratic kernels the agammaegation is stronger when the contribution from the diffusion term is small. However, similarly to the polynomial kernel a connected mass is not maintained in either of these cases. Additionally, in the cases where a connected mass is maintained we still observe a sharp drop off in the middle of the profile - indicating formation of a separation. We believe that the tendency to separate in this case is caused by a combination of the strong local interaction and the strong contribution from the food sources; particles that start close to a food source are "trapped" and remain close to the food source as the contribution from the food source and the local interaction reinforce each other. Then, particles towards the center of the domain feel a pull towards food sources form both the food source itself and the "trapped" particles.

4.2.   Interaction dominated regime

We now examine simulations of Eq (2.17) in the regime where the interaction term is the strongest contribution. Again, we start with the case where the kernel function $W$ is quadratic and given by $W(x) = \frac{x^{2}}{2}$. The results can be seen in Figure 6. Here, we do not see any food source agammaegation. This can possibly be explained by the fact that in these simulations the initial profile is Gaussian. We have previously seen that Gaussians are stationary states for Eq (2.17) in the case of a quadratic kernel and no drift term.

Next we examine the case where the kernel function is polynomial, of the same form as in the food source dominated regime. The results can be seen in Figure 7. Here, analogous to the food source dominated regime we see more pronounced agammaegation in the cases where the diffusion term has the weakest contribution. However the agammaegation is less pronounced than the respective examples in the food source dominated regime due to the stronger contribution from the interaction term. Unlike the analogous cases in the drift dominated regime we do not observe food source agammaegation in the cases where $C = 5$, this is possibly due to the stronger contribution from the short range repulsion in the interaction term.

Finally, we examine the evolution in the case of a Gaussian kernel function. The results can be seen in Figure 8. Here, in contrast to the respective cases in the drift dominated regime, we do not observe food source agammaegation in the cases where the diffusion term has the second largest contribution. Addtionally, unlike the analogous cases in the drift dominated regime and the polynomial kernel in the interaction dominated regime, we do not observe agammaegation in the case where $B = 10, C = 1, A = 1$. We do observe agammaegation in the case where $A = 5, B = 10, C = 1$. However, similarly to the polynomial kernel the effect is less pronounced than in the analogous case in the drift dominated regime. Both of these effects are consistent with a stronger contribution from the attractive kernel.

4.3.   Competition regime

We now simulate cases where two parameters dominate in competition. The results can be seen in Figures 9-11. Interestingly, across all three interaction kernels, pronounced food source agammaegation is only observed in cases where $A = B = 10$. All of these cases except one exhibit the qualitative features we wish to capture as they do maintain a connected mass. All other cases do not result in food source agammaegation. This can be explained by the fact that in all of these cases the diffusion term has a significant contribution. The observation that strong contributions from the diffusion term dampen food source agregation is a general pattern that we have also observed in the other regimes and is additionally supported by the statement of Theorem 2.

The cases in the competition regime where we do see food source agammaegation also reinforce another pattern observed across all regimes; examples that have the qualitative features we wish to capture have strong contributions from the drift term, the interaction term, or both. However, we also find that cases in which a connected mass was not maintained only occured where the drift term was dominant. Likewise, cases in which food source agammaegation did not occur in the absence of strong diffusion only occured in the interaction dominated regime. These observations suggest a general modeling principal; contributions from the drift and diffusion terms should be scaled similarly and larger than the contribution from the diffusion term.

4.4.   Evidence for asymptotic states and a two dimensional example

We now turn to some further numerical studies aimed at validating observations made in the previous section. First, note that in the previous section we used a uniform simulation time. Therefore, in theory some of the qualitative observations we have made could merely be artifacts of not simulating the dynamics for long enough; its possible that we are not observing asymptotic states of the dynamics. For example, in Figure 3 we observe what appears to be food source agammaegation with a connected mass in the regime $A = 10, B = 5, C = 1$ and mass separation in the regime $A = 10, B = 1, C = 1$ suggesting that the stronger interaction term in the first case helps the profile maintain connectivity. However, if we are not truly observing asymptotic states its possible that we would observe mass separation if we waited a longer time. In Figure 12 we perform the same simulation as Figure 3 except the simulation time is doubled and we plot six profiles instead of four. However, most of the six profiles are not distinguishable as they are overlapping in exactly the same profiles as the final profiles plotted in Figure 3. Mass agammaegation with a connected profile is observed in exactly analogous fashions in the regime $A = 10, B = 5, C = 1$ and mass separation is observed in both cases in the regime $A = 10, B = 1, C = 1$. This provides strong evidence that we are observing asymptotic states, however to provide further support in this direction (and to show that this phenomena does not depend on the kernel choice) we repeat all simulations of the polynomial kernel with the longer simulation time in Figures 13-15. We again observe exactly analogous results to the shorter simulation time in all cases providing further support that asymptotic states are observed.

The aim of this study was to identify parameter regimes in which the derived macroscopic model qualitatively displays slime mold food seeking behavior. For this reason we have limited our simulations to one dimension as they are faster to compute and if a phenomena does not appear in one dimension it will surely not manifest in two. However, for actual application to slime mold movement (and to possibly validate the model with empirical experiments) two dimensions is surely the case of interest. To illustrate that observations made in one dimension do manifest in two dimensions we will show that mass separation observed in one dimension is also observed in two dimensions as we feel that this phenomena is easiest to qualitatively evaluate in the two dimensional plots. We simulate the macroscopic model in two dimensions with the gaussian kernel in the regime $A = 10, B = 5, C = 1$ - the results can be found in Figure 16. In this case the dynamics were simulated in the region $[-2.1, 2.1]\times[-2.1, 2.1]$, the food sources are located at the points $(1, 0)$ and $(-1, 0)$, 225 particles were used in the blob method ODE, and the initial profile is given by the two dimensional gaussian centered at the origin with the same variance as that used in one dimension. Analogous to the one dimensional case we observe mass separation.

5.   Conclusions

In this manuscript we presented a model hierarchy aimed at modeling the food seeking behavior of the slime mold P. polycephalum. We were principally interested in capturing the slime mold's ability to agammaegate around disparate food sources while maintaining a connected mass. We first presented a particle based model which includes three main features - a drift term to model a gradient of chemoattractant produced by food sources, an interaction term to model the slime mold's propensity to maintain a connected mass, and a diffusion term to model slime mold foraging behavior. Additionally, we included scaling parameters for each of these three terms in order to study how varying their contributions could aid in modeling P. polycephalum. However simulating enough particles to realistically model the evolution of slime mold would be computationally intractable. Therefore, under the assumption of propagation of chaos we showed that in the large particle limit the evolution of the particle model can be described by a macroscopic agammaegation-diffusion equation which can be efficiently simulated. Before embarking on simulating the equation we made some analytical observations via analysis of the equation's stationary states. First, we found that in the case of no drift term that an assumption of Gaussian stationary states implies that the interaction kernel must be quadratic and discussed how this could be used as a heuristic to determine a kernel that realistically models P. polycephalum behavior. Then, we showed that if the contribution from the diffusion term is sufficiently larger than the interaction and drift terms that the only possible stationary state is 0, deeming this regime as unsuitable for modeling P. polycephalum.

We then ran a series of simulations of the macroscopic model in order to investigate how to scale the three contributing terms and choose the interaction kernel in order to model P. polycephalum food seeking behavior. We simulated three different parameter regimes; a regime in which the drift term dominated, a regime in which the interaction term dominated and a regime in which the dominating terms were in competition. In each regime we examined three different interaction kernels, an attractive kernel with a local interaction, an attractive kernel with long range interaction, and a kernel that included a short range repulsion. We observed cases that reproduced connected mass food source agammaegation across all three regimes with all kernels however several general patterns did emerge. We found that strong contributions from the diffusion term tend to dampen food source agammaegation even in the presence of strong contributions from the other terms. Additionally, we found that in the absence of strong diffusion, dominance of the drift term can result in mass separation while dominance of the interaction term can prevent food agammaegation. This suggests a general modeling strategy of scaling contributions from the drift and interaction terms similarly and larger than the diffusion term.

There are many ways in which this study could be extended. Agammaegation-diffusion equations have been widely studied in the case of no drift term ^[27]. A main finding is that there is a critical mass above which the interaction term becomes dominant and causes a "blow up" in the mass profile. From the particle view this can be thought of intuitively as all particles converging on one point. Below the critical mass, the tendency of the diffusion term to cause the mass profile to spread prevents a blow up from occuring. We did not observe this phenomena in any of our simulations (even in the interaction dominated regime). It would be interesting to consider whether the presence of a drift term can prevent blow up from occuring at all or mitigate it in the sense that the critical mass at which blow up occurs becomes higher. Additionally, as we did observe the desired qualitative behavior in several cases it would be useful to compare the evolution of our model to data from actual slime mold growth in a simple situation such as the two food source case we consider. This would be an important first step in the quantitative validation of this model.

Acknowledgments

Part of this work was carried out while D. Weber was visiting the University of Mannheim within the program IPID4all funded by the German Academic Exchange Service (DAAD). This work was supported by the DAAD project Stochastic dynamics for complex networks and systems (Project-ID 57444394).

The publication of this article was funded by the Baden-Württemberg Ministry of Science, Research and the Arts and by the University of Mannheim.

Conflict of interests

All authors declare no conflicts of interest in this paper.