Mathematical Modelling of Turning Delays in Swarm Robotics

We investigate the effect of turning delays on the behaviour of groups of differential wheeled robots and show that the group-level behaviour can be described by a transport equation with a suitably incorporated delay. The results of our mathematical analysis are supported by numerical simulations and experiments with e-puck robots. The experimental quantity we compare to our revised model is the mean time for robots to find the target area in an unknown environment. The transport equation with delay better predicts the mean time to find the target than the standard transport equation without delay.


Introduction
Much theory has been developed for the coordination and control of distributed autonomous agents, where collections of robots are acting in environments in which only short-range communication is possible (Reif and Wang, 1999).By performing actions based on the presence or absence of signals, algorithms have been created to achieve some greater group level task; for instance, to reconnoitre an area of interest whilst collecting data or maintaining formations (Desai et al., 2001).In this paper, we will investigate an implementation of searching algorithms, similar to those used by flagellated bacteria, in a robotic system.
Many flagellated bacteria such as Escherichia coli (E.coli) use a run-and-tumble searching strategy in which movement consists of more-or-less straight runs interrupted by brief tumbles (Berg, 1983).When their motors rotate counter-clockwise the flagella form a bundle that propels the cell forward with a roughly constant speed; when one or more flagellar motors rotate clockwise the bundle flies apart and the cell 'tumbles' (Kim et al., 2003).Tumbles reorient the cell in a more-or-less uniformly-random direction (with a slight bias in the direction of the previous run) for the next run (Berg and Brown, 1972).In the absence of signal gradients the random walk is unbiased, with a mean run time ∼ 1 sec and a tumble time ∼ 0.1 sec.However, when exposed to an external signal gradient, the cell responds by increasing (decreasing) the run length when moving towards (away from) a favorable direction, and therefore the random walk is biased with a drift in that direction (Berg, 1975;Koshland, 1980).Similar behaviour can be observed in swarms of animals avoiding predators and coordinating themselves within a group (Couzin et al., 2002).
The behaviour of E. coli is often modelled as a velocity-jump process where the time spent tumbling is neglected as it is much smaller than the time spent running (Othmer et al., 1988;Erban and Othmer, 2004).In such a velocity-jump process, particles follow a given velocity u from a set of allowed velocities V ⊂ R d , d = 2, 3, for a finite time.The particle changes velocity probabilistically according to a Poisson process with intensity λ , i.e. the mean run-duration is 1/λ .A new velocity v is chosen Assuming that λ and T are constant, one can show that the long-time behaviour of the density ̺(t, x) = V p(t, x, v) dv is given by the diffusion equation (Hillen and Othmer, 2000).If λ depends on an external signal (e.g.nutrient concentration), then the resulting velocity jump process is biased and its long time behaviour can be described by a drift-diffusion equation for ̺ (Othmer and Hillen, 2002;Erban and Othmer, 2005).
In this paper, we will study an experimental system based on E-Puck robots (Bonani and Mondada, 2004).We programme these differential wheeled robots to follow a run-and-tumble searching strategy in order to find a given target set.We will concentrate on the simplest possible scenario: an unbiased velocity-jump process in two spatial dimensions with the fixed speed s ∈ R + , the constant mean run time λ −1 ∈ R + , and the turning kernel which is independent of u A special feature of the E-Puck robots is that they can perform turns on the spot as in the classical velocity-jump process described by (1.1).In this paper, we will investigate in how far (1.1) presents a good description of the behaviour of the robotic system and we will develop an extension of (1.1) that results in a better match between experimental data and mathematical model.The paper is organized as follows: in Section 2, we introduce the experimental system as well as the obtained data.This data is compared to the classical velocity-jump theory.In Section 3, we extend the velocity-jump theory to include finite turning times and we compare it to our experimental data, showing a much improved match.We conclude our paper, in Section 4, by discussing the implications of our results .

Velocity-jump processes in experiments with robots
Equation (1.1) introduced the density behaviour of the general velocity-jump process that we are aiming to investigate using the experimental set-up described in Section 2.1.In particular, we will concentrate on a simple unbiased velocity-jump process with the fixed speed s ∈ R + , the mean run duration λ −1 ∈ R + and the turning kernel (1.2).Whilst throughout this work, we keep λ constant through time, it is common to model λ in biological applications as varying in response to external signals or internal variables (Erban and Othmer, 2005).This fixed-speed unbiased velocity-jump process can be viewed as a starting point for considering more complex searching algorithms.We will demonstrate that by including a small modification (the introduction of a delay to the turning kernel), we can alter this simple velocity-jump process so that it models the behaviour of the E-Puck robots.
We are interested in comparing the idealised velocity-jump process, given in (1.1)-(1.2) to robotic experiments.Due to a restriction in numbers of robots, one cannot feasibly talk about a "density" of robots that could be compared to p(t, x, v) as given in (1.1).Therefore, our experiments concentrate on the escape of robots from a given domain.We may interpret this as the target finding ability of the E-Puck robots.Using these experiments, we can infer data both on the flux at the barrier and the exit times and can compare those to numerical results of velocity-jump processes in Sections 2.3 and 2.4.

Experimental set-up and procedure
To obtain the empirical data, an experimental system consisting of 16 E-Puck robots was used.E-Puck robots are small differential wheeled robots with a programmable microchip (Bonani and Mondada, 2004).The diameter of each robot is ε = 75 mm with a height of 50 mm and weight of 200 g.Throughout the experiments, the speed was chosen to be s = 5.8 × 10 −2 m/sec.The robots turn with an angular velocity ω = 4.65/sec.Full specifications along with a picture are given in Appendix A.
In the experiments, we use a rectangular arena Ω with walls on three of the 4 edges and an opening to the target area T along the fourth edge1 .A diagram of the arena along with the notation used can be seen in Figure 1 and a photo is shown in Figure 5(b) in Appendix A. When considering such an arena, one has to distinguish between the size of the physical arena and the effective arena (shown in blue in Figure 1) that the robot centres can occupy.The effective arena used in the experiments has the dimensions L x = 1.183 m and L y = 1.145 m = L x − ε/2.The reflective (wall) boundary and the target boundary will be denoted as ∂ Ω R and ∂ Ω T , respectively, and can be defined as Throughout the remainder of the paper, we will also use n R (resp.n T ) to denote the outwards pointing normal on the reflective (resp.target) boundary.
During the experiments, robots were initialised inside a removable square pen Ω 0 of effective edge length L 0 = 0.305 m, shown in Figure 1 and Figure 5(b) in Appendix A. A short period of free movement within the pen before its removal allowed us to reliably release all robots into the full domain Ω at the same time as well as randomising their initial positions within the pen.We recorded the exit time for each of the robots, when its geometric centre entered the target area T .Each repetition of the experiment was continued for 300 sec or until all 16 robots had left the arena.
The robots were programmed using C and a cross-compiling tool, with the firmware being transferred onto the robots via bluetooth.A pseudo-code of the algorithm implemented on the robots is shown in Table 1.This algorithm represents a velocity-jump process in the limit as ∆t → 0 (Erban et al., 2006), and gives a good approximation as long as λ ∆t ≪ 1.In the experiments we used λ = 0.25 sec −1 implying a mean run duration of 4 sec and ∆t = 0.1 sec, resulting in λ ∆t = 2.5 × 10 −2 ≪ 1.Note that, s, ω and λ can be changed on a software level on an E-Puck.For ω we chose the maximum possible value, whilst for s we chose a value below the physical maximum.Choosing a lower velocity means that we mitigate the effects of acceleration and deceleration to the running speed since the robots cannot do this instantaneously as the basic velocity jump model assumes.In a practical setting, one could interpret s and ω as given characteristics of the system, whilst λ can be chosen in a way that accelerates the target finding process for the given application with the choice of λ likely to represent a trade-off between sampling an area and time spent reorienting.
In addition to the algorithm in Table 1, robots were also made to implement an obstacle-avoidance strategy the using the four proximity sensors placed at angles ±17.5 • and ±47.5 • from the centre axis in the front part of the E-Puck.Reflective turns were carried out based on the signals received at these sensors.As the robots are incapable of distinguishing between walls and other robots, those reflections occur whether a robot collides with the wall ∂ Ω R or another robot.As a consequence we discuss the importance of robot-robot collisions on the experimental results in Section 2.2.

Relevance of collisions for low numbers of robots
For non-interacting particles which can change direction instantaneously, equation (1.1) accurately describes the mesoscopic density through time.However, in our experiments the robots undergo reflective collisions when they come into close contact, rather than passing through or over each other.For a low number of particles, we used Monte Carlo simulations to demonstrate that collisions are not the dominant behaviour and have little effect on the distribution of particles.In Figure 2, we compare two Monte Carlo simulations: (a) in which particles are allowed to pass through one another and (b) in which collisions are modelled explicitly.In (c) we present the solution of equation (1.1).This comparison demonstrates that the mean density of the underlying process converges to the solution of the integro-differential equation (IDE).The parameters employed in this model comparison are taken directly from the equivalent robot experiment; (s, λ , ε) = (5.8× 10 −2 m/sec, 0.250 sec −1 , 7.5 × 10 −2 m).In Figure 2(c), for the differential equation, we use a first-order numerical scheme with ∆ θ = π/20, ∆ x = L y /200, ∆t = 10 −2 sec.
In the Monte Carlo simulations we initialise particles in the effective pen for 20 sec where they undergo hard-sphere collisions.They are then released into the larger arena where in one simulation they are point-particles and in the other they undergo reflective collisions as hard-spheres.Instead of [1] Robot is started at position x(0) ∈ Ω 0 .Generate r 1 ∈ [0, 1] uniformly at random, set t = 0 and v(0) = s cos 2πr 1 sin 2πr 1 .
Table 1.An algorithmic implementation of the velocity jump process.
(a) removing particles at the target boundary as shown in Figure 1 (as we do in the robot experiments), this edge of the domain is closed so that all edges correspond to reflective boundary conditions.For the IDE (1.1), we model the initial condition as a step function over the pen.These densities are visualised in Figure 2. Formally, this initial condition can be written as where χ Ω 0 denotes the indicator function of the initial region Ω 0 .The corresponding boundary condition is p(t, x, v) = p(t, x, v ′ ) for x ∈ ∂ Ω R where the reflected velocity v ′ is defined as where n Ω is the outward pointing normal at the position x ∈ ∂ Ω .
After 20 sec, we record the density in each of the scenarios and present the results in Figure 2.There is minimal visible discrepancy between the Monte Carlo simulations presented in Figure 2 for our choice of parameter values.In order to compare the three simulations given in Figure 2 we also employed a pairwise Kolmogorov-Smirnov test (Peacock, 1983).A value (of the Kolmogorov-Smirnov metric) close to zero denotes a good fit between the two simulations.It corresponds to the probability that one can reject the hypothesis that the distributions are identical.When comparing the two Monte Carlo simulations, a value of 2.37 × 10 −2 was obtained; when comparing equation (1.1) with the hard-sphere Monte Carlo simulation, a value of 5.65 × 10 −2 was obtained; finally when comparing equation (1.1) with the point-particle Monte Carlo simulation, a value of 3.40 × 10 −2 was obtained.This demonstrates that all three density distributions are all highly similar.
In the limit where N → ∞, for N being the number of hard-spheres, the IDE (1.1) can be altered by the addition of a Boltzmann-like collision term (Harrisi, 1971;Cercignani, 1988).We will report various density limits of hard-sphere velocity-jump processes in a future publication.

Comparison between theory and experiments: loss of mass over time
In this and subsequent sections we compare the results of 50 repetitions of the experiments described in Section 2.1 with numerical results obtained by solving the corresponding mathematical equations.One way of interpreting the experimental exit-time data is by considering the expected mass remaining inside the arena Ω at a given time.For the experimental data this quantity is plotted as a solid line in Figure 3(a).We compare this result to the variation of the remaining mass with time from a numerical solution of (1.1) combined with the following boundary conditions: where the reflected velocity v ′ is defined by (2.3).As demonstrated in Section 2.2, such a comparison is reasonable since collisions do not have a major impact in the parameter regime chosen here.The initial condition for the PDE (1.1) is identical to the condition given in equation (2.2).The mass remaining in the domain in the PDE system is then defined as and is plotted as a dashed line in Figure 3(a).The initial mass is normalized to 1.
As an additional verification of the validity of (1.1) for this system as well as the negligibility of robot-robot interactions discussed in Section 2.2, we performed Monte Carlo simulations with and without collisions and both gave a result almost indistinguishable from the dashed line in Figure 3(a).An obvious observation from Figure 3(a) is that the PDE description does not match the experimental data well, with the robots exiting the arena significantly slower than predicted.In this figure, we use a firstorder finite volume method with ∆ θ = π/20, ∆ x = 1.183 m/200 and ∆t = 10 −3 sec in order to solve the IDE (1.1).

Comparison between theory and experiments: mean exit time problem
An alternative way to interpret the experimental data is to consider mean exit times.We calculate the mean time taken for a robot to leave the arena and plot this value as the solid line in Figure 3(b).In order to be able to compare this value to analytic results, one has to reformulate the transport equation (1.1) into a mean exit time problem.Consider the conditional mean exit time for a particle that it leaves the arena before t end , i.e. the expected time until exit, given that it leaves before the experiment ends.To determine the mean exit time, τ = τ(x 0 , v 0 ), we specify starting position x 0 and velocity v 0 .As we are using conditional mean exit time problems, we introduce ρ end , which is the probability that the particle is still in Ω at the end time t end .The mean exit time satisfies the following equation (2.5) In Section 3, in which delays are modelled, a derivation is given for the mean exit time problem; setting the delay term to zero allows one to see how equation (2.5) is derived.This so-called "backwards problem" satisfies the following boundary conditions where v ′ 0 is again the reflected velocity with respect to v 0 as defined in (2.3).Due to the arena shape, by taking the spatial average in the y-direction one can further simplify the mean exit time problem.In the case where the turning kernel is given by equation (1.2), one can obtain a problem with two parameters x 0 and θ , where θ ∈ (−π, π] is the angle defining the velocity v 0 by v 0 = s(cos(θ ), sin(θ )).For τ x = τ x (x 0 , θ ) (2.8) Of the 50 × 16 = 800 robots in the experiments, 92 did not leave the arena before time t end = 300 sec.Therefore, we use the value ρ end ≈ 0.115.When initial direction cannot be specified, the mean-exit time from a given x-position is given by 1 2π where τ x is the solution of (2.8).This is plotted as the dashed line in Figure 3(b).The numerical solution was performed using an upwind-scheme in the x-direction with ∆ x = 1.1825 m/200 and an angular discretisation of ∆ θ = π/20.Additionally, we take the spatial average of the mean-exit time from the initial region Ω 0 and plot this as the bold dashed line in Figure 3(b).This line does not match well with the corresponding average mean-exit time found in the robot experiments.The mean exit time of the robots was 137.1 sec, whilst the numerical solution of equation (2.8) predicted 120.1 sec, meaning an underestimation of 17.0 sec or 14.3%.In the following section we will extend the classical velocity-jump theory to improve this match with the experimental data.

Modelling turning delays
In Section 2.2, we observed that collisions between robots does not play a major role in explaining the discrepancy between the transport equation (1.1) and the experimental data presented in Sections 2.3 and 2.4.As well as assuming independently moving particles, the transport equation (1.1) is also predicated on the assumption that the reorientation phase takes a negligible amount of time compared to the running phase.Since this assumption is not satisfied in our robot experiments, this section extends the original model through the inclusion of finite turning times.

Introduction of a resting state
Let us initially state two assumptions that apply to the robot experiment, but might not extend to velocity jump processes in biological systems, like the run-and-tumble motion of E. Coli (Berg, 1983), which has motivated the searching strategies implemented on robots: (a) a new direction v ′ ∈ V is chosen as soon as the particle enters the reorientation ("tumble") phase; (b) the time it takes for a particle to reorient ("tumble") from velocity Assumption (b) implies that the turning time is constant in time and equal for each particle and, in particular, does not depend on the particle's history.For the robots studied in this paper, we can additionally assume that reorientation phase is equivalent to a directed rotation with a constant angular velocity ω ∈ R + .Therefore, the turning time depends only on the angle between the current velocity v ∈ V and the new velocity v ′ ∈ V and K takes the form We now extend the classical model (1.1) through the introduction of a resting state r(t, x, v, η) that formally defines the number of particles currently "tumbling" (turning) towards their new chosen velocity v and remaining turning time η.The density p(t, x, v) will now only denote the particles which are at time t in the run phase.The update of the extended system is given through In (3.2) we can see that running particles will enter a tumble phase with rate λ and particles that have finished the tumble signified through η = 0 will re-enter the run-phase.Equation (3.3) represents the linear relation between t and η and shows that particles enter the tumble phase depending on their newly chosen velocity direction.In order to guarantee conservation of mass throughout the system, we introduce the non-negativity condition for η through Additionally, the boundary conditions for the system (3.2)-(3.3) are given through where v ′ is the reflected velocity of v given by (2.3).In order to show that the system (3.2)-(3.3) is actually consistent, we prove that mass in the system is conserved if no target is present.
LEMMA 3.1 The total mass in system (3.2)-(3.3)with the boundary conditions given in (3.4) in the case of reflective boundaries everywhere Proof.We define for every point x ∈ ∂ Ω R the two subsets V + and V − of V as follows Additionally, let us define Integrating (3.3) with respect to η ∈ [0, ∞), we obtain after reordering for x / ∈ ∂ Ω : Hence, for every point x / ∈ ∂ Ω we obtain Integrating this with respect to x ∈ Ω and v ∈ V gives Using the divergence theorem, we can evaluate the integral on the right hand side to be where we have used the second boundary condition in (3.4) in the last step.Additionally, for x ∈ ∂ Ω and v ∈ V − (x), we obtain by integrating the third boundary condition in (3.4) with respect to η Integrating this with respect to x ∈ ∂ Ω and v ∈ V and using the last boundary condition in (3.4) we obtain Summing up the results from (3.6) and (3.7), we obtain dM/dt = 0 and hence the total mass M(t) in the system is conserved.

Transport equation with turning delays
We eliminate the resting state from system (3.2)-(3.3)and derive the generalization of the transport equation (1.1) to a transport equation with a suitably incorporated delay.This can be done by solving (3.3) for r using the method of characteristics, which results in where H is the Heaviside step function.Let us assume that K(v, u) is given by (3.1).Then K(v, u) π/ω.Considering times t > π/ω, we have r(0, x, v,t) = 0. We can now substitute (3.8) into (3.2) to for t > π/ω.Note that (3.9)only considers particles in the running phase and hence does not strictly conserve mass.The boundary conditions for transport equation (3.9) are (3.10)

Equation for mean-exit time
Equation (3.9) can be rewritten as M p = 0, where the operator M is given by For a forward problem specified by M p = 0, coupled with initial and boundary conditions, the backward problem is given by the adjoint operator M * q = 0, with final condition and adjoint boundary conditions (Lewins, 1965).The adjoint operator is given by: Using integration by parts and the divergence theorem, we see where we used the boundary conditions We will also assume the following boundary conditions (3.14) Then the last term in (3.12) is equal to zero as it is shown in Appendix B. Using (3.12)-(3.14)and the variable set (t 0 , x 0 , v 0 ) to indicate starting times and positions, we can write the backwards equation M * q = 0 in the following form: (3.15)More precisly, we should write q(t 0 , x 0 , v 0 ) ≡ p(t, x, v |t 0 , x 0 , v 0 ), i.e. q gives the probability that the particle is at the position x with velocity v at time t given that its initial position and velocity at time t 0 were x 0 and v 0 , respectively.Let ρ ≡ ρ(t, x 0 , v 0 ) be the probability that the particle is in Ω at time t given that the initial position and velocity is given as x 0 and v 0 , respectively.Then Substituting t 0 = −t into (3.15) and using the Taylor expansion, we obtain ∂ ρ ∂t The probability of a single particle leaving where we denoted the probability that a particle is in Ω at the end time by ρ end = ρ(t end , x 0 , v 0 ).Integrating (3.16) over time, we obtain (3.17) where we neglected the higher order terms.By Taylor-expanding the boundary terms from equation (3.10) and integrating in time, we obtain the following boundary conditions where the reflected velocity v ′ 0 is given by (2.3), i.e.

Comparison between the transport equation theory with delays and experimental results
Let us now compare the extended theory developed in Sections 3.1-3.3 to the experimental data using the same approach as in Sections 2.3 and 2.4.For the case of the arena given in Figure 1, we write Ω = (0, L x ) × (−L y /2, L y /2) and T = (L x , ∞) × (−L y /2, L y /2) and we simplify equation (3.17) by integrating over the y-direction to obtain an average value for τ for our position along the x-axis.Let us define this average: 0 ), integrating (3.17) and using (3.18), we obtain the following equation for τ (3.20) In the case where T is the unbiased, fixed-speed, 2-dimensional turning kernel given by (1.2) and using (3.1), we have v 0 = (v (x) 0 , v (y) 0 ) = s(cos θ , sin θ ) and we can evaluate the second integral term in equation (3.20) explicitly to be Then (3.20) can be rewritten as follows where A(θ ) is defined by The boundary conditions (3.18) simplify to The numerical solution of (3.21)-(3.22)can be further simplified by considering the symmetry in angle τ x (x 0 , θ ) = τ x (x 0 , −θ ), i.e. it is sufficient to solve (3.21)where (x 0 , θ ) are restricted to the domain (0, L x ) × (0, π) with boundary conditions (3.22).

Comparison between theory and experiments: loss of mass over time
In this section, we show that the transport theory with delays better explains the experimental data with robots by considering the loss of mass over time, as we did in Section 2.3.In Figure 4(a), for the same experimental set-up as in Figure 3, we plot the mass remaining in the system against time.The solid line represents the experimental data, whilst the classical theory from Figure 3(a) is re-plotted using the dashed line.The red dash-dotted line (color online) now shows a numerical solution of system (3.2)-(3.3)that incorporates the finite reorientation time into the analysis.The numerical solution was achieved using a first order finite-volume method paired with an upwind scheme for (3.3).For (3.2) we used ∆ x = 1.183 m/200, ∆t = 10 −3 sec and ∆ θ = π/20.For (3.3) we used the same ∆t = 10 −3 sec and a discretisation of ∆ η = 3.38 × 10 −2 sec corresponding to the time it takes to turn from one velocity direction to the next.Figure 4(a) demonstrates that the inclusion of turning delays provides an improved match to the experimental data.

Comparison between theory and experiments: mean exit time problem
The mean exit time problem from Section 2.4 can also be better modelled by the transport equation theory with suitable incorporated delays as it is demonstrated in Figure 3 This numerical solution was obtained using the same method as in Section 2.4 and we again plot the average over the initial pen as a bold dashed line.The bold dashed line indicates a predicted mean exit time of 133.2 sec compared to the experimental value 137.3 sec, an error of approximately 3.1% or 4.1 sec.This represents a strong improvement to the discrepancy of 14.3% seen for the model that neglected the turning events (dashed line) and goes to show that turning times are indeed non-negligible and can be built into our model in a consistent manner.

Discussion
In this paper, we have studied an implementation of a run-and-tumble searching strategy in a robotic system.The algorithm implemented by the robots is motivated by a biological system -behaviour of the flagellated bacterium E. coli.Bio-inspired algorithms are relatively common in swarm robotics.Algorithms based on behaviour of social insects have been implemented previously in the literature, see for example Garnier et al. (2005); Krieger et al. (2000); Webb (2000) and Fong et al. (2003).One of the challenges of bio-inspired algorithms is that robots do not have the same sensors as animals.
For example, E. coli bias their movement according to extracellular chemical systems.In biological models, chemical signals often evolve according to the solution of reaction-diffusion partial differential equations (Franz and Erban, 2012;Franz et al., 2013).Therefore, an implementation of the full run-andtumble chemotactic model in the robotic system requires either special sensors for detecting chemical signals, e.g.robots for odor detecting (Russell, 2001), or replacing chemical signals by suitable caricatures of them, e.g. using glowing floor for E-Puck robots (Mayet et al., 2010).
The main goal of this paper is to study how the mathematical theory developed for E. coli applies to the robotic system based on E-Pucks.Thus we do not focus on technological challenges connected with sensing chemical signals or their analogues (Russell, 2001;Mayet et al., 2010) and the implemented velocity jump process is unbiased, given by the turning kernel (1.2).In this case the mathematical theory is well understood.If the collisions between particles (robots or bacteria) and reorientation times can be neglected, then this velocity jump process is described by the transport equation (1.1) and the long time behaviour is given by the diffusion equation (Hillen and Othmer, 2000).In Section 2.2 we show that collisions between robots are negligible in our experimental set up.However, we still observe quantitative differences between the results based on the transport equation (1.1) and robotic experiments.
In Section 3 we identify turning delays as the main mechanism contributing to differences between the mathematical theory developed for E. coli and the results of experiments with E-Pucks.We introduce the resting state in equations (3.2)-(3.3)and then derive the transport equation with delay (3.9).Our delay term is different from models of tumbling of E. coli, because the underlying physical process is different.Tumbling times of E. coli are exponentially distributed, i.e. they can be explicitly included in mathematical models by using transport equations which take into account probabilistic changes to and from the resting (tumbling) state (Erban and Othmer, 2004).In the case of robots, the turning time depends linearly on the turning angle.The selection of new direction is effectively instant and the main contributing factor to turning delays is the finite turning speed of robots.
We have studied a relatively simple searching algorithm motivated by E. Coli behaviour, but the transport equations and velocity jump processes naturally appear in modelling of other biological systems, such as modelling chemotaxis of amoeboid cells (Erban and Othmer, 2007) or swarming behaviour as seen in various fish, birds and insects (Carrillo et al., 2009;Erban and Haskovec, 2012).We conclude that the same delay terms as in (3.9) would be applicable whenever we implement these models in E-Pucks.From a mathematical point of view, it is also interesting to consider coupling of (3.9) with sensing of extracellular signals, because signal transduction also has its own delay which can be modelled using velocity jump models with internal dynamics (Franz et al., 2013;Erban and Othmer, 2004;Xue and Othmer, 2009).We are also currently investigating how the transport equations for velocityjump processes need to be adapted for higher densities of robots.We will report our findings in a future publication.

FIG. 1 .
FIG. 1. Schematic showing of the experimental set-up along with the notation used throughout this paper.Dotted border lines correspond to the effective arena and bold lines to the actual arena.For further details see the text.
FIG. 2. Comparison of individual-based simulations with (1.1).Each plot shows the resulting density at the final time of the simulation, 20 sec.(a) Individual-based simulation using point particles.(b) Individual-based simulation using hard-sphere interactions.(c) Numerical solution to (1.1) using a finite volume method with parameters given in the text.
FIG. 3. (Colour available online.)Comparing the experimental results to the solution of the IDE (1.1).Panel (a) displays the evolution of the mean mass remaining in the system.The solid line (black) corresponds to the experimental data and the dashed line (blue) to the numerical solution of (1.1) with boundary conditions (2.4).Panel (b) displays the mean exit time averaged over all velocities vs the initial distance from the absorbing boundary.The solid line (black) represents the experimental data and the dashed line (blue) the mean exit time derived from the numerical solution of (1.1) with boundary conditions (2.6).In order to allow direct comparison with the experimental data, the shorter bold dashed line represents the average of theoretically derived exit time (thin dashed line) over the region Ω 0 , from which the robots were released in the experimental scenario.For both plots parameters and numerical methods are described in the text.
FIG. 4. (Colour available online.)Comparing the experimental results (solid black line) to the solution of the system of differential equations (3.2)-(3.3).Panel (a) as Figure 3(a), the dash-dotted line (red) shows the numerical solution of the system of equations (3.2)-(3.3)with boundary conditions (3.4).Panel (b) as Figure 3(b), the dash-dotted (red) line shows the mean exit time computed using equation (3.21) with boundary conditions (3.22).In order to allow direct comparison with the experimental data, the shorter bold lines represents the average of theoretically derived exit times over the region Ω 0 , from which the robots were released in the experimental scenario.The previous plots from Figure3are shown as the dashed line (blue).For both plots parameters and numerical methods are described in the text.
(b).The solid line represents again the experimental data, whilst the classical results are shown as dashed lines.The numerical solution of (3.21) with the boundary conditions (3.22) is shown as the red dot-dashed line (color online).