File preview
Modeling Populations of Micro-robots for Biological Applications
Paul Bogdan, Guopeng Wei, and Radu Marculescu
Department of Electrical and Computer Engineering Carnegie Mellon University, Pittsburgh, PA, USA {pbogdan, guopengw, radum}@ece.cmu.edu
Abstract— In order to detect and cure diseases affecting microscopic regions of the human body, swarms of micro-robots can harness the motility of bacteria and swim towards inaccessible regions of the body to detect abnormal behavior and/or deliver various drugs. In this paper, we propose a modeling approach for the dynamics of dense networks of swarms and estimate their performance via stochastic metrics. Keywords-component; Random walks, collective dynamics, competitive behavior, dense networks of swarms, disease diagnosis, drug delivery.
swimming in the spinal cord and sensing the environment for detecting potential cancer risks. II. RELATED WORK AND NOVEL CONTRIBUTION Understanding and modeling of biological propelled microrobotic swarms has its roots in biological studies of cooperation and competition between various organisms such as bacteria, fish, birds, or herds of mammals. Generally speaking, the mathematical models of biological swarming can be classified in two classes: individual-based (aka Langevin or Langrangian) [1][6][7][9][18][19][22] and continuum-based models [2][8] [16][18][24]. Individual-based modeling relies on motion equations for each individual of the swarm. In contrast, continuum-based modeling describes the swarm dynamics via advection-diffusion reaction equations where the swarm state variable is the concentration or population density. In addition to these models, several swarm modeling approaches deal with partial differential equation characterizing generic (non-space dependent) metrics, such as the distribution of group sizes or motion orientation [10][21][25]. In spite of this significant body of work related to biological and micro-robotic swarming, there are still important problems that need to be solved. One such problem is related to the modeling of micro-robotic swarms, while taking into account the characteristics of both long- and short-range microscopic interactions between swarm agents. This mathematical challenge is critical not only for obtaining more accurate macroscopic dynamic models, but also to better understand related phenomena at cellular level. A secondary but equally important problem is to construct non-Fickian models which avoid the limitations of point-like models and account for the spatial volume occupied by micro-robots. To address the above mentioned challenges and to study the dynamics of dense networks of micro-robotic swarms, we represent micro-robots dynamics as a collection of interacting random walks which obey the volume exclusion rules [13] in a 3D space. The new 3D multiple random walk formalism takes into-account agent-to-agent interactions and sets forth the basis for accurate description of microscopic swarm dynamics. III. MODELING OF DENSE NETWORKS OF SWARMS In order to characterize the dynamics of micro-robotic swarms, we map the movement of multiple micro-robots within a confined 3D region to the problem of tracking the trajectories of multiple simultaneous random walks which can interact at
O
I.
NE
INTRODUCTION
of the main challenges of modern medicine is the detection of silent progression and spreading of various diseases through the human body. For instance, cancer is one of the leading causes of death because in most situations, tumors appear and develop undetected by many of the screening tests or even by the immune system itself. Indeed, while a large body of research suggests that tumors silent progression is frequently accompanied by an increased demand of nutrients and oxygen [12], these events are not easily detected by the immune system and the cancer cells can corrupt the neighboring and remote organs up to the point where surgery cannot help with a cure. Despite these challenges, we can rely on bacteria motility [3] to engineer micro-robotic swarms capable of monitoring, detection and targeted drug delivery within the human body [20]. However, while fabrication of single micro-robot is well understood, the mathematical characterization of the collective dynamics of swarms [9][15][19] of such micro-robots still represents a major challenge from the perspective of modeling, analyzing and designing such dense swarms of bacteria for targeted diagnostic and drug delivery.
In this paper, we propose a mathematical model for capturing the dynamics of a large number (or teams) of selfdriven micro-robots (i.e., bacteria propelled capsules) able to swim and access small regions of the human body; this voyage takes place in a non-invasive manner due to micro-robots dimensions [20]. Such engineered micro-robots can perform massively parallel and distributed tasks related to diagnostic or drug delivery purposes. Given the affinity of chemotactic bacteria (e.g., Serratia Marcescens) to high oxygen consumption around tumors, the micro-robots can sense and swim through the interstitial spaces towards the affected regions. For instance, Fig. 1 shows micro-robotic swarms
Fig. 1. Dense network of micro-robotic swarms swim in hard to access regions of the body. The micro-robots dynamics is modeled as a collective behavior of multiple interacting random walks in a 3D graph tessellation of space. Each graph node has associated a binary random variable (σ) denoting whether or not it is occupied by one micro-robot. For detection and health monitoring purposes the goal is to find the timedependent coverage of the swarm. The goal of the drug delivery problem is to find the probability for the nodes in the disease area (see Fig.1.b) to be occupied by micro-robots.
a)
Micro-robot at location h on the 3D tessellation graph
b) L
Targeted drug delivery region: Find the probability of σm(t) = 1, ∀ m ∈Critical_Area
σ p-1qw = 1
σ pqw = 1
No micro-robot at location i-1jk
σ i-1jk = 0
Starting configuration
Σij+1k = 1 σijk = 1
Final configuration
various points in space and time. Towards this end, we first tessellate the 3D space into a graph of size M×N×Q and, as shown in Fig. 1.a, associate with each node (i,j,k) a binary random variable σijk which indicates whether or not there exists a micro-robot at location (i,j,k) (i.e., σijk=1 if the micro-robot is present and σijk=0 if it’s not). In order to capture the realm of biological world, we also assume that the graph tessellation of space is fine-grained (i.e., in the order of a micro-robot size) and impose exclusion conditions such that a graph node can be occupied by only a single micro-robot. In addition, we consider the chemical communication that can exist between bacteria and model the interactions between neighboring micro-robots. To study the dynamics of multiple interacting random walks contained within a certain region, we need to determine the probability P(σ111,…, σijk,…,σMNQ;t) of having multiple random walkers at arbitrary locations in the 3D lattice. The time evolution of probability P(σ111,…, σijk,…,σMNQ;t) is dictated by the following transition probabilities: a) Moving backwards along X-axis: We assume that, during an infinitesimal time interval δt, a random walker (i.e., micro-robot) located at (i,j,k) location is moving to the (i-1,j,k) location and satisfies the following transition probability:
Pr σ i −1 jk (t + δt ) + 1,σ ijk (t + δt ) | σ i −1 jk (t ),σ ijk (t ) + 1 = α ijk →i −1 jk σ ijk (t ) + 1 1 − σ i −1 jk (t ) δt + O(δt )
where σijk(t) and σi+1jk(t) are binary random variables representing whether locations (i,j,k) and (i+1,j,k) are occupied by a random walker at time t, α ijk → i +1 jk is the transition rate. c) Moving backwards along Y-axis: During a short time interval δt, a random walker at (i,j,k) location on a graph can move to the (i,j-1,k) location according to the following transition probability:
= α ijk →ij −1k σ ijk (t ) + 1 1 − σ ij −1k (t ) δt + O(δt ) Pr σ ij −1k (t + δt ) + 1,σ ijk (t + δt ) | σ ij −1k (t ),σ ijk (t ) + 1
{
(
)(
)
}
(3)
where σijk(t) and σij-1k(t) are binary random variables representing whether locations (i,j,k) and (i,j-1,k) are occupied by a random walker at time t, α ijk → ij −1k is the transition rate. d) Moving forward along Y-axis: The random walker at (i,j,k) location on a graph can move to the (i+1,j,k) location according to the following transition probability:
= α ijk →ij +1k σ ijk (t ) + 1 1 − σ ij +1k (t ) δt + O(δt ) Pr σ ijk (t + δt ),σ ij +1k (t + δt ) + 1 | σ ijk (t ) + 1,σ ij +1k (t )
{
(
)(
)
}
(4)
{
(
)(
)
}
(1)
where σijk(t) and σij+1k(t) are binary random variables representing whether locations (i,j,k) and (i,j+1,k) are occupied by a random walker at time t, α ijk → ij +1k is the transition rate. e) Moving backwards along Z-axis: Similarly, in a short time interval δt, a random walker at (i,j,k) location on a graph can move to the (i,j,k-1) location according to the following transition probability:
Pr σ ijk −1 (t + δt ) + 1,σ ijk (t + δt ) | σ ijk −1 (t ),σ ijk (t ) + 1 = α ijk →ijk −1 σ ijk (t ) + 1 1 − σ ijk −1 (t ) δt + O(δt )
where σijk(t) and σi-1jk(t) are binary random variables representing whether locations (i,j,k) and (i-1,j,k) are occupied by a random walker (micro-robot) at time t, α ijk → i −1 jk is the transition rate and O(δt) shows that all terms higher than δt can be neglected. To better understand the intuition behind (1), let us consider that σi-1jk(t)=1 at time t. In this case, because the location (i-1,j,k) is occupied by a random walker, the transition equation in (1) is zero and so the random walker at (i,j,k) cannot move to (i-1,j,k) location. b) Moving forward along X-axis: During an infinitesimal time interval δt, a random walker (micro-robot) at any (i,j,k) location on a graph can move to the (i+1,j,k) location according to the following transition probability:
= α ijk →i +1 jk σ ijk (t ) + 1 1 − σ i +1 jk (t ) δt + O(δt ) Pr σ ijk (t + δt ),σ i +1 jk (t + δt ) + 1 | σ ijk (t ) + 1,σ i +1 jk (t )
{
(
)(
)
}
(5)
where σijk(t) and σijk-1(t) are binary random variables representing whether locations (i,j,k) and (i,j,k-1) are occupied by a random walker at time t, α ijk →ijk −1 is the transition rate. f) Moving forward along Z-axis: In a short time interval δt, a random walker at (i,j,k) location on a graph can move to the (i,j,k+1) location according to the following transition probability:
{
(
)(
)
}
(2)
Fig. 2. Random walkers (micro-robots) in a 3D space. A detailed description of such micro-robots is given in [3] a) Represent the 3D space as a lattice tessellation with a lattice edge of 1 µm. Assign to each node (i,j,k) of the lattice a time-dependent binary random variable σijk(t) which becomes 1 (or 0) when a micro-robot is present (or absent) at this node or around a sphere of 1µm radius centered at this particular node. b) We assume the micro-robots are injected 1 cm away from the targeted region. These micro-robots propelled by chemotactic bacteria sense the target and swim to the target following a bias random walk.
a)
b)
σijk+1 = 0 σi-1jk = 0 σij-1k = 0 σijk = 1 σi+1jk = 0 σijk-1 = 1 σij+1k = 0
Targeted (destination) locations for the microrobots
Insertion (source) locations for the microrobots
Pr σ ijk (t + δt ),σ ijk +1 (t + δt ) + 1 | σ ijk (t ) + 1,σ ijk +1 (t ) = α ijk →ijk +1 σ ijk (t ) + 1 1 − σ ijk +1 (t ) δt + O(δt )
{
(
)(
)
}
(6)
where σijk(t) and σijk+1(t) are binary random variables representing whether locations (i,j,k) and (i,j,k+1) are occupied by a random walker at time t, α ijk →ijk +1 is the transition rate. Taking into account these transition probabilities, the evolution of the probability P(σ111,…, σijk,…,σMNQ;t) is given by the following master equation:
dP σ111,...,σ ijk ,...,σ MNQ ; t dt
M , N ,Q i , j , k =1
For tumor detection and monitoring purposes, the focus of this analysis is on finding the coverage and frequency of visiting certain nodes in the 3Ds space by at least one random walker. This can be done by appropriately choosing a metric that averages over a certain space region and specific time interval the joint probability P(σ111,…, σijk,…,σMNQ;t) determined by (7). In contrast, solving the drug delivery problem requires finding the probability of having a critical number of microrobots within a specific area (see the highlighted region in Fig. 2.b). Note that classical diffusion theory [17] cannot be applied to such a scenario since the dynamics of such multiple random walks are crucially affected by various hydro-dynamical and chemical interactions. This non-equilibrium statistical physics approach is meant to capture the collective and competitive behavior of particles and predict the evolution of the swarm as a function of the density of walkers and the strength of their interactions. As shown, the accurate modeling of trajectories and distances travelled by micro-robotic swarms is of crucial importance for solving the diagnostic and drug delivery problems. IV. SIMULATION RESULTS Solving (7) in 3D is not a trivial task, since it implies to determine the state probabilities in a space consisting of a large number of states. Consequently, in this section, we simulate the master equation in (7) and record the trajectories of T interacting random walks to investigate their implications on the drug delivery problem. To get more intuition about the mathematical form of the probability P(σ111,…,σijk,…,σMNQ;t), we simulate 1000 interacting random walkers starting from 1000 random locations on the lower left hand side (as shown in Fig. 2.b) and moving only along the edges towards 1000 random targets on the upper right hand side of a 3D 104×104×104 lattice which corresponds to a space of 104×104×104 μm. .In the absence of neighboring micro-robots (i.e., no interaction), any micro-robot has an equal probability in choosing any of the six possible directions with a small bias of 0.016 toward target sources. This small bias is meant to model the affinity of micro-robots to move towards regions exhibiting an abnormal behavior and which can create chemical gradients. Generally speaking, this bias should be related to the amount of abnormal chemical gradients and the response mechanism of micro-robots to neighboring stimulus.
(
)=
∑ {α
i −1 jk →ijk
(1 − σ ijk )σ i −1 jk P(σ111,...,σ ijk ,...,σ MNQ ; t ) +
( ) ( ) ( ) ( ) (7) + α ij +1k →ijk (1 − σ ijk )σ ij +1k P (σ111,...,σ ijk ,...,σ MNQ ; t ) + + α ijk −1→ijk (1 − σ ijk )σ ijk −1P (σ 111,...,σ ijk ,...,σ MNQ ; t ) + + α ijk +1→ijk (1 − σ ijk )σ ijk +1P (σ 111,...,σ ijk ,...,σ MNQ ; t ) − σ ijk ⋅ [αijk →i −1 jk (1 − σ i −1 jk ) + αijk →i +1 jk (1 − σ i +1 jk ) +αijk →ij −1k ⋅ (1 − σ ij −1k ) + αijk →ij +1k (1 − σ ij +1k ) + αijk →ijk −1(1 − σ ijk −1 ) + α ijk →ijk +1 (1 − σ ijk +1 )]P (σ111,...,σ ijk ,...,σ MNQ ; t )}
+ α i +1 jk →ijk 1 − σ ijk σ i +1 jk P σ111,...,σ ijk ,...,σ MNQ ; t + + α ij −1k →ijk 1 − σ ijk σ ij −1k P σ111,...,σ ijk ,...,σ MNQ ; t +
For completeness, the probability P(σ111,…, σijk,…,σMNQ;t) must satisfy the following conservation laws:
∑∑∑σ ∑ P(σ
i =1 j =1 k =1 Q
ijk = 0,1
M
N
Q
111,..., σ ijk ,..., σ MNQ ; t
)=1
(8)
∑∑∑σ ∑σ P(σ
ijk i =1 j =1 k =1
ijk = 0,1
M
N
111,..., σ ijk ,..., σ MNQ ; t
) = NRW
(9)
where NRW in eq.(9) represents the number of random walkers (micro-robots) under investigation. Of note, although the above mathematical formalism has been written based on a one-hop interaction between random walkers, it can be generalized for arbitrary distances by re-defining the transition probabilities in (1) to (6) to account for specific interaction radius. This is left for future work.
Fig. 3. a ) Probability of visiting a sphere (denoted by a node with 3D coordinates) of 5 µm radius in a 104×104×104 volume by any of the biased micro-robots at time t exhibits an almost bell-shaped behavior. Higher interaction distance among micro-robots causes the probability of visiting a certain sphere to be more spread over time. b) Number of distinct spheres covered by at multiple random walkers as a function of the interaction radii in a 104×104×104 volume.
x 10
6
Probability of visiting a node at time t
Number of distinct spheres covered p
a)
Node at (50,50,50) 1 hop interaction Gaussian fitting Node at (50,50,50) 15 hop interaction Gaussian fitting
8 7 6 5 4 3 2 1 0 2
Non-interacting RWs non-interacting RWs 20 hops interaction RWs - 20 hops interaction 25 hops interaction RWs - 25 hops interaction 10 hops interaction RWs - 10 hops interaction
0
b)
Simulation time [abstract units]
Simulation time [abstract units] Simulation time [abstract units]
4
6
8
10
12
14
16
x 10
4
In addition, when two or more micro-robots are within a predefined interaction radius, their motion is biased away from each other and so they tend to cover more space and increase the likelihood of detecting abnormalities. These interactions can be implemented by means of molecular communication mechanisms as described in [11][14][23]. Consequently, although the initial bias forces the micro-robots to drift toward a few cohesive swarms, the predefined interactions made them to repel each other for small intervals of time. To ensure statistical significance, we extract the results presented herein from 1000 independent simulations of micro-robotic populations that took up to 7 days of computational runtime on a 16 core Intel Xeon computer with 16GB of physical memory. Fig. 3a shows that the probability of visiting a certain location in the 3D 104×104×104 space by any of the interacting micro-robots (as a function of time) exhibits a sparse bell shaped behavior. By increasing the interaction radius from 1 hop to 15 hops causes the probability of visiting a certain location to be more spread over time. This implies that interaction radius is a very important parameter in enhancing the probability of detecting any abnormal behavior by increasing the likelihood that a certain location is visited over extended periods of time. By considering interaction radius of 10, 20 and 25 hops between neighboring micro-robots, we measure the efficiency of the dense networks of micro-robotic swarms in detecting cancer cells via the number of distinct spheres covered by time t. The results shown in Fig. 3b were obtained by averaging over 1000 simulation runs. Fig. 3b shows that using the power of multiple random walkers and exploiting the interaction radius between microrobots can contribute to a significant coverage in a short interval of time. However, due to spatial constraints this is not conserved at medium and long times. From Fig. 3b one can notice that the swarms exhibit a non-monotonic behavior in terms of number of distinct covered/visited spheres as a function of the interaction radius. This has significant implications on the design of interacting micro-robotic swarms. Indeed, by designing the interacting micro-robots independently of the spatial constraints can lead to a fast search in the early times but also to a saturation point in terms of the region visited by the micro-robots and this can lead to inefficient detection of any abnormal behavior. For the current experimental setup, one can safely assume that an interaction
radius of 20 hops offers the best detection of the abnormal behavior. Nevertheless, in practical cases the interaction radius need to be carefully set as a function of the space in which the micro-robots move while considering the characteristics of the movement process (i.e., the directional transition probabilities). Besides coverage, it is also essential to quantify how much time it takes for a group of micro-robots to reach a certain targeted region; this would be a hitting time metric [5]. Consequently, we measured the hitting times for all microrobots to reach their predefined destinations for four interaction radius: 1 hop (1 μm), 5 hops (5 μm), 10 hops (10 μm), 15 hops (15 μm), respectively. Fig. 4 shows the probability of hitting time to exceed a certain threshold as a function the four interaction radius. The hitting time probabilities were also obtained from 1000 independent simulation runs. One can observe that by considering larger interaction radii, the probability of hitting time to exceed a given threshold tends diverge from the Gaussian law and is better approximated by long tail type of distributions such as Fisk, Student’s t and generalized extreme value. V. CONCLUSIONS AND FUTURE WORK In this paper, we have set forth a non-equilibrium statistical physics description of micro-robotic swarms with application to abnormal behavior detection and drug delivery. Unlike previous swarm models, we have developed a mathematical formalism which takes into consideration the volume exclusion rules and micro-robot-to-micro-robot interactions. To get more intuition about the implications of spatial dimensions of microrobots and their interactions, we have presented a comprehensive simulation investigation. Our results show that volume exclusion rules and interaction radii have significant implications on the probability of detecting abnormal behavior (e.g., cancer cells, fungal infections) at both short and long times. In addition, volume exclusion rules and interaction radii also affect the distribution of drug delivery (hitting) times. Consequently, when designing dense networks of microrobotic swarms it is essential to take into account the dimensions of micro-robots, their volume exclusion rule, and the nature of their interactions. ACKNOWLEDGMENT This work was supported by the US National Science Foundation (NSF) under Grant CPS-1135850.
a)
b)
c)
d)
Fig. 4 Probability (Hitting times> threshold) for various interaction distances among micro-robots within a 1cm×1cm×1cm region. a) 1 hop interaction. b) 5 hops interaction. c) 10 hops interaction. d) 15 hops interaction. In scenarios c) and d), the distribution of hitting times is better approximated by means of the generalized extreme value, Fisk and Student’s t-distributions rather than by means of the Gaussian law.
References
R. Axelrod, The Complexity of Cooperation: Agent-Based Models of Competition and Collaboration, Princeton University Press, 1997. [2] R. E. Baker, C. A. Yates and R. Erban, “From Microscopic to Macroscopic Descriptions of Cell Migration on Growing Domains,” Bulletin of Mathematical Biology, Vol. 72, No. 3, pp. 719-762, 2010. [3] B. Behkam and M. Sitti, “Characterization of Bacterial Actuation of Micro-Objects,” IEEE Intl. Conf. on Robotics and Automation, 2009. [4] H. C. Berg, Random Walks in Biology, Princeton Univ. Press, 1993. [5] P. Bogdan and R. Marculescu, “Hitting Time Analysis for Fault-Tolerant Communication at Nanoscale in Future Multiprocessor Platforms,” IEEE Trans. on CAD of Integrated Circuits and Systems, vol.30, no.8, 2011 [6] C. M. Breder, “Equations Descriptive of Fish Schools and Other Animal Aggregations,” Ecology, Vol. 35, No. 3, pp. 361-370, Jul., 1954. [7] D. L. Deangelis and L. J. Gross (Editors), Individual-Based Models and Approaches in Ecology: Populations, Communities and Ecosystems, Chapman and Hall. 1992. [8] A. Galstyan, T. Hogg, and K. Lerman, “Modeling and Mathematical Analysis of Swarms of Microscopic Robots,” IEEE Swarm Intelligence Symposium (SIS-2005), Pasadena, CA, 2005. [9] V. Gazi and K.M. Passino, Swarm Stability and Optimization, Springer Technology & Engineering, 2011. [10] S. Gueron and S. A. Levin, “The Dynamics of Group Formation,” Mathematical Biosciences, Vol. 128, Issues 1–2, pp. 243-264, 1995. [11] E. Gul, B. Atakan, and O.B. Akan, “NanoNS: A Nanoscale Network Simulator Framework for Molecular Communications”, Nano Communication Networks Journal (Elsevier), vol. 1, pp. 138-156, 2010. [12] D. Hanahan and R.A. Weinberg, “The Hallmarks of Cancer,” Cell, Vol. 100, Issue 1, pp. 57-70, January 2000. [1]
[13] C. Kipnis and C. Landim, Scaling Limits of Interacting Particle Systems, Springer - Mathematics, 1999. [14] I. Llatser, I. Pascual, N. Garralda, A. Cabellos-Aparicio, and E. Alarcón, “N3Sim: A Simulation Framework for Diffusion-based Molecular Communication”, IEEE TC on Simulation, no. 8, pp. 3-4, March 2011. [15] R. Mahnke, J. Kaupužs and I. Lubashevsky, Physics of Stochastic Processes: How Randomness Acts in Time, Wiley Mathematics, 2009. [16] A.Mogilner and L.Edelstein-Keshet, “A Non-Local Model for a Swarm,” Journal of Mathematical Biology, vol. 38 , pp. 534-570, 1999. [17] J. D. Murray, Mathematical Biology, Springer - Science, 2003. [18] G. Naldi, Mathematical Modeling of Collective Behavior in SocioEconomic and Life Sciences, Springer-Science, 2010. [19] F. Schweitzer and J. D. Farmer, Brownian Agents and Active Particles: Collective Dynamics in the Natural and Social Sciences, Springer, 2007. [20] M. Sitti, “Miniature Devices: Voyage of the microrobots,” Nature, vol. 458, pp. 1121-1122, 30 April 2009. [21] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, “Novel Type of Phase Transition in a System of Self-Driven Particles,” Phys. Rev. Lett. , vol. 75, pp. 1226–1229, 1995. [22] K. Warburton and J. Lazarus, “Tendency-Distance Models of Social Cohesion in Animal Groups,” Journal of Theoretical Biology, 1991. [23] R. Weiss and T.F. Knight, Jr. “Engineered Communications for Microbial Robotics,” 6th Intl. Workshop on DNA-Based Computers: DNA Computing, A. Condon, G. Rozenberg (Eds.) Springer, 1-16, 2000. [24] C. Xue and H. G. Othmer, “Multiscale Models of Taxis-Driven Patterning in Bacterial Populations,” SIAM Journal on Applied Mathematics, vol. 70, pp. 133-167, 2009. [25] C.A. Yates, R.E. Baker, R. Erban, and P.K. Maini, “Refining SelfPropelled Particle Models for Collective Behaviour,” Canadian Journal of Applied Mathematics, vol. 18, pp. 299-350, 2010.
Paul Bogdan, Guopeng Wei, and Radu Marculescu
Department of Electrical and Computer Engineering Carnegie Mellon University, Pittsburgh, PA, USA {pbogdan, guopengw, radum}@ece.cmu.edu
Abstract— In order to detect and cure diseases affecting microscopic regions of the human body, swarms of micro-robots can harness the motility of bacteria and swim towards inaccessible regions of the body to detect abnormal behavior and/or deliver various drugs. In this paper, we propose a modeling approach for the dynamics of dense networks of swarms and estimate their performance via stochastic metrics. Keywords-component; Random walks, collective dynamics, competitive behavior, dense networks of swarms, disease diagnosis, drug delivery.
swimming in the spinal cord and sensing the environment for detecting potential cancer risks. II. RELATED WORK AND NOVEL CONTRIBUTION Understanding and modeling of biological propelled microrobotic swarms has its roots in biological studies of cooperation and competition between various organisms such as bacteria, fish, birds, or herds of mammals. Generally speaking, the mathematical models of biological swarming can be classified in two classes: individual-based (aka Langevin or Langrangian) [1][6][7][9][18][19][22] and continuum-based models [2][8] [16][18][24]. Individual-based modeling relies on motion equations for each individual of the swarm. In contrast, continuum-based modeling describes the swarm dynamics via advection-diffusion reaction equations where the swarm state variable is the concentration or population density. In addition to these models, several swarm modeling approaches deal with partial differential equation characterizing generic (non-space dependent) metrics, such as the distribution of group sizes or motion orientation [10][21][25]. In spite of this significant body of work related to biological and micro-robotic swarming, there are still important problems that need to be solved. One such problem is related to the modeling of micro-robotic swarms, while taking into account the characteristics of both long- and short-range microscopic interactions between swarm agents. This mathematical challenge is critical not only for obtaining more accurate macroscopic dynamic models, but also to better understand related phenomena at cellular level. A secondary but equally important problem is to construct non-Fickian models which avoid the limitations of point-like models and account for the spatial volume occupied by micro-robots. To address the above mentioned challenges and to study the dynamics of dense networks of micro-robotic swarms, we represent micro-robots dynamics as a collection of interacting random walks which obey the volume exclusion rules [13] in a 3D space. The new 3D multiple random walk formalism takes into-account agent-to-agent interactions and sets forth the basis for accurate description of microscopic swarm dynamics. III. MODELING OF DENSE NETWORKS OF SWARMS In order to characterize the dynamics of micro-robotic swarms, we map the movement of multiple micro-robots within a confined 3D region to the problem of tracking the trajectories of multiple simultaneous random walks which can interact at
O
I.
NE
INTRODUCTION
of the main challenges of modern medicine is the detection of silent progression and spreading of various diseases through the human body. For instance, cancer is one of the leading causes of death because in most situations, tumors appear and develop undetected by many of the screening tests or even by the immune system itself. Indeed, while a large body of research suggests that tumors silent progression is frequently accompanied by an increased demand of nutrients and oxygen [12], these events are not easily detected by the immune system and the cancer cells can corrupt the neighboring and remote organs up to the point where surgery cannot help with a cure. Despite these challenges, we can rely on bacteria motility [3] to engineer micro-robotic swarms capable of monitoring, detection and targeted drug delivery within the human body [20]. However, while fabrication of single micro-robot is well understood, the mathematical characterization of the collective dynamics of swarms [9][15][19] of such micro-robots still represents a major challenge from the perspective of modeling, analyzing and designing such dense swarms of bacteria for targeted diagnostic and drug delivery.
In this paper, we propose a mathematical model for capturing the dynamics of a large number (or teams) of selfdriven micro-robots (i.e., bacteria propelled capsules) able to swim and access small regions of the human body; this voyage takes place in a non-invasive manner due to micro-robots dimensions [20]. Such engineered micro-robots can perform massively parallel and distributed tasks related to diagnostic or drug delivery purposes. Given the affinity of chemotactic bacteria (e.g., Serratia Marcescens) to high oxygen consumption around tumors, the micro-robots can sense and swim through the interstitial spaces towards the affected regions. For instance, Fig. 1 shows micro-robotic swarms
Fig. 1. Dense network of micro-robotic swarms swim in hard to access regions of the body. The micro-robots dynamics is modeled as a collective behavior of multiple interacting random walks in a 3D graph tessellation of space. Each graph node has associated a binary random variable (σ) denoting whether or not it is occupied by one micro-robot. For detection and health monitoring purposes the goal is to find the timedependent coverage of the swarm. The goal of the drug delivery problem is to find the probability for the nodes in the disease area (see Fig.1.b) to be occupied by micro-robots.
a)
Micro-robot at location h on the 3D tessellation graph
b) L
Targeted drug delivery region: Find the probability of σm(t) = 1, ∀ m ∈Critical_Area
σ p-1qw = 1
σ pqw = 1
No micro-robot at location i-1jk
σ i-1jk = 0
Starting configuration
Σij+1k = 1 σijk = 1
Final configuration
various points in space and time. Towards this end, we first tessellate the 3D space into a graph of size M×N×Q and, as shown in Fig. 1.a, associate with each node (i,j,k) a binary random variable σijk which indicates whether or not there exists a micro-robot at location (i,j,k) (i.e., σijk=1 if the micro-robot is present and σijk=0 if it’s not). In order to capture the realm of biological world, we also assume that the graph tessellation of space is fine-grained (i.e., in the order of a micro-robot size) and impose exclusion conditions such that a graph node can be occupied by only a single micro-robot. In addition, we consider the chemical communication that can exist between bacteria and model the interactions between neighboring micro-robots. To study the dynamics of multiple interacting random walks contained within a certain region, we need to determine the probability P(σ111,…, σijk,…,σMNQ;t) of having multiple random walkers at arbitrary locations in the 3D lattice. The time evolution of probability P(σ111,…, σijk,…,σMNQ;t) is dictated by the following transition probabilities: a) Moving backwards along X-axis: We assume that, during an infinitesimal time interval δt, a random walker (i.e., micro-robot) located at (i,j,k) location is moving to the (i-1,j,k) location and satisfies the following transition probability:
Pr σ i −1 jk (t + δt ) + 1,σ ijk (t + δt ) | σ i −1 jk (t ),σ ijk (t ) + 1 = α ijk →i −1 jk σ ijk (t ) + 1 1 − σ i −1 jk (t ) δt + O(δt )
where σijk(t) and σi+1jk(t) are binary random variables representing whether locations (i,j,k) and (i+1,j,k) are occupied by a random walker at time t, α ijk → i +1 jk is the transition rate. c) Moving backwards along Y-axis: During a short time interval δt, a random walker at (i,j,k) location on a graph can move to the (i,j-1,k) location according to the following transition probability:
= α ijk →ij −1k σ ijk (t ) + 1 1 − σ ij −1k (t ) δt + O(δt ) Pr σ ij −1k (t + δt ) + 1,σ ijk (t + δt ) | σ ij −1k (t ),σ ijk (t ) + 1
{
(
)(
)
}
(3)
where σijk(t) and σij-1k(t) are binary random variables representing whether locations (i,j,k) and (i,j-1,k) are occupied by a random walker at time t, α ijk → ij −1k is the transition rate. d) Moving forward along Y-axis: The random walker at (i,j,k) location on a graph can move to the (i+1,j,k) location according to the following transition probability:
= α ijk →ij +1k σ ijk (t ) + 1 1 − σ ij +1k (t ) δt + O(δt ) Pr σ ijk (t + δt ),σ ij +1k (t + δt ) + 1 | σ ijk (t ) + 1,σ ij +1k (t )
{
(
)(
)
}
(4)
{
(
)(
)
}
(1)
where σijk(t) and σij+1k(t) are binary random variables representing whether locations (i,j,k) and (i,j+1,k) are occupied by a random walker at time t, α ijk → ij +1k is the transition rate. e) Moving backwards along Z-axis: Similarly, in a short time interval δt, a random walker at (i,j,k) location on a graph can move to the (i,j,k-1) location according to the following transition probability:
Pr σ ijk −1 (t + δt ) + 1,σ ijk (t + δt ) | σ ijk −1 (t ),σ ijk (t ) + 1 = α ijk →ijk −1 σ ijk (t ) + 1 1 − σ ijk −1 (t ) δt + O(δt )
where σijk(t) and σi-1jk(t) are binary random variables representing whether locations (i,j,k) and (i-1,j,k) are occupied by a random walker (micro-robot) at time t, α ijk → i −1 jk is the transition rate and O(δt) shows that all terms higher than δt can be neglected. To better understand the intuition behind (1), let us consider that σi-1jk(t)=1 at time t. In this case, because the location (i-1,j,k) is occupied by a random walker, the transition equation in (1) is zero and so the random walker at (i,j,k) cannot move to (i-1,j,k) location. b) Moving forward along X-axis: During an infinitesimal time interval δt, a random walker (micro-robot) at any (i,j,k) location on a graph can move to the (i+1,j,k) location according to the following transition probability:
= α ijk →i +1 jk σ ijk (t ) + 1 1 − σ i +1 jk (t ) δt + O(δt ) Pr σ ijk (t + δt ),σ i +1 jk (t + δt ) + 1 | σ ijk (t ) + 1,σ i +1 jk (t )
{
(
)(
)
}
(5)
where σijk(t) and σijk-1(t) are binary random variables representing whether locations (i,j,k) and (i,j,k-1) are occupied by a random walker at time t, α ijk →ijk −1 is the transition rate. f) Moving forward along Z-axis: In a short time interval δt, a random walker at (i,j,k) location on a graph can move to the (i,j,k+1) location according to the following transition probability:
{
(
)(
)
}
(2)
Fig. 2. Random walkers (micro-robots) in a 3D space. A detailed description of such micro-robots is given in [3] a) Represent the 3D space as a lattice tessellation with a lattice edge of 1 µm. Assign to each node (i,j,k) of the lattice a time-dependent binary random variable σijk(t) which becomes 1 (or 0) when a micro-robot is present (or absent) at this node or around a sphere of 1µm radius centered at this particular node. b) We assume the micro-robots are injected 1 cm away from the targeted region. These micro-robots propelled by chemotactic bacteria sense the target and swim to the target following a bias random walk.
a)
b)
σijk+1 = 0 σi-1jk = 0 σij-1k = 0 σijk = 1 σi+1jk = 0 σijk-1 = 1 σij+1k = 0
Targeted (destination) locations for the microrobots
Insertion (source) locations for the microrobots
Pr σ ijk (t + δt ),σ ijk +1 (t + δt ) + 1 | σ ijk (t ) + 1,σ ijk +1 (t ) = α ijk →ijk +1 σ ijk (t ) + 1 1 − σ ijk +1 (t ) δt + O(δt )
{
(
)(
)
}
(6)
where σijk(t) and σijk+1(t) are binary random variables representing whether locations (i,j,k) and (i,j,k+1) are occupied by a random walker at time t, α ijk →ijk +1 is the transition rate. Taking into account these transition probabilities, the evolution of the probability P(σ111,…, σijk,…,σMNQ;t) is given by the following master equation:
dP σ111,...,σ ijk ,...,σ MNQ ; t dt
M , N ,Q i , j , k =1
For tumor detection and monitoring purposes, the focus of this analysis is on finding the coverage and frequency of visiting certain nodes in the 3Ds space by at least one random walker. This can be done by appropriately choosing a metric that averages over a certain space region and specific time interval the joint probability P(σ111,…, σijk,…,σMNQ;t) determined by (7). In contrast, solving the drug delivery problem requires finding the probability of having a critical number of microrobots within a specific area (see the highlighted region in Fig. 2.b). Note that classical diffusion theory [17] cannot be applied to such a scenario since the dynamics of such multiple random walks are crucially affected by various hydro-dynamical and chemical interactions. This non-equilibrium statistical physics approach is meant to capture the collective and competitive behavior of particles and predict the evolution of the swarm as a function of the density of walkers and the strength of their interactions. As shown, the accurate modeling of trajectories and distances travelled by micro-robotic swarms is of crucial importance for solving the diagnostic and drug delivery problems. IV. SIMULATION RESULTS Solving (7) in 3D is not a trivial task, since it implies to determine the state probabilities in a space consisting of a large number of states. Consequently, in this section, we simulate the master equation in (7) and record the trajectories of T interacting random walks to investigate their implications on the drug delivery problem. To get more intuition about the mathematical form of the probability P(σ111,…,σijk,…,σMNQ;t), we simulate 1000 interacting random walkers starting from 1000 random locations on the lower left hand side (as shown in Fig. 2.b) and moving only along the edges towards 1000 random targets on the upper right hand side of a 3D 104×104×104 lattice which corresponds to a space of 104×104×104 μm. .In the absence of neighboring micro-robots (i.e., no interaction), any micro-robot has an equal probability in choosing any of the six possible directions with a small bias of 0.016 toward target sources. This small bias is meant to model the affinity of micro-robots to move towards regions exhibiting an abnormal behavior and which can create chemical gradients. Generally speaking, this bias should be related to the amount of abnormal chemical gradients and the response mechanism of micro-robots to neighboring stimulus.
(
)=
∑ {α
i −1 jk →ijk
(1 − σ ijk )σ i −1 jk P(σ111,...,σ ijk ,...,σ MNQ ; t ) +
( ) ( ) ( ) ( ) (7) + α ij +1k →ijk (1 − σ ijk )σ ij +1k P (σ111,...,σ ijk ,...,σ MNQ ; t ) + + α ijk −1→ijk (1 − σ ijk )σ ijk −1P (σ 111,...,σ ijk ,...,σ MNQ ; t ) + + α ijk +1→ijk (1 − σ ijk )σ ijk +1P (σ 111,...,σ ijk ,...,σ MNQ ; t ) − σ ijk ⋅ [αijk →i −1 jk (1 − σ i −1 jk ) + αijk →i +1 jk (1 − σ i +1 jk ) +αijk →ij −1k ⋅ (1 − σ ij −1k ) + αijk →ij +1k (1 − σ ij +1k ) + αijk →ijk −1(1 − σ ijk −1 ) + α ijk →ijk +1 (1 − σ ijk +1 )]P (σ111,...,σ ijk ,...,σ MNQ ; t )}
+ α i +1 jk →ijk 1 − σ ijk σ i +1 jk P σ111,...,σ ijk ,...,σ MNQ ; t + + α ij −1k →ijk 1 − σ ijk σ ij −1k P σ111,...,σ ijk ,...,σ MNQ ; t +
For completeness, the probability P(σ111,…, σijk,…,σMNQ;t) must satisfy the following conservation laws:
∑∑∑σ ∑ P(σ
i =1 j =1 k =1 Q
ijk = 0,1
M
N
Q
111,..., σ ijk ,..., σ MNQ ; t
)=1
(8)
∑∑∑σ ∑σ P(σ
ijk i =1 j =1 k =1
ijk = 0,1
M
N
111,..., σ ijk ,..., σ MNQ ; t
) = NRW
(9)
where NRW in eq.(9) represents the number of random walkers (micro-robots) under investigation. Of note, although the above mathematical formalism has been written based on a one-hop interaction between random walkers, it can be generalized for arbitrary distances by re-defining the transition probabilities in (1) to (6) to account for specific interaction radius. This is left for future work.
Fig. 3. a ) Probability of visiting a sphere (denoted by a node with 3D coordinates) of 5 µm radius in a 104×104×104 volume by any of the biased micro-robots at time t exhibits an almost bell-shaped behavior. Higher interaction distance among micro-robots causes the probability of visiting a certain sphere to be more spread over time. b) Number of distinct spheres covered by at multiple random walkers as a function of the interaction radii in a 104×104×104 volume.
x 10
6
Probability of visiting a node at time t
Number of distinct spheres covered p
a)
Node at (50,50,50) 1 hop interaction Gaussian fitting Node at (50,50,50) 15 hop interaction Gaussian fitting
8 7 6 5 4 3 2 1 0 2
Non-interacting RWs non-interacting RWs 20 hops interaction RWs - 20 hops interaction 25 hops interaction RWs - 25 hops interaction 10 hops interaction RWs - 10 hops interaction
0
b)
Simulation time [abstract units]
Simulation time [abstract units] Simulation time [abstract units]
4
6
8
10
12
14
16
x 10
4
In addition, when two or more micro-robots are within a predefined interaction radius, their motion is biased away from each other and so they tend to cover more space and increase the likelihood of detecting abnormalities. These interactions can be implemented by means of molecular communication mechanisms as described in [11][14][23]. Consequently, although the initial bias forces the micro-robots to drift toward a few cohesive swarms, the predefined interactions made them to repel each other for small intervals of time. To ensure statistical significance, we extract the results presented herein from 1000 independent simulations of micro-robotic populations that took up to 7 days of computational runtime on a 16 core Intel Xeon computer with 16GB of physical memory. Fig. 3a shows that the probability of visiting a certain location in the 3D 104×104×104 space by any of the interacting micro-robots (as a function of time) exhibits a sparse bell shaped behavior. By increasing the interaction radius from 1 hop to 15 hops causes the probability of visiting a certain location to be more spread over time. This implies that interaction radius is a very important parameter in enhancing the probability of detecting any abnormal behavior by increasing the likelihood that a certain location is visited over extended periods of time. By considering interaction radius of 10, 20 and 25 hops between neighboring micro-robots, we measure the efficiency of the dense networks of micro-robotic swarms in detecting cancer cells via the number of distinct spheres covered by time t. The results shown in Fig. 3b were obtained by averaging over 1000 simulation runs. Fig. 3b shows that using the power of multiple random walkers and exploiting the interaction radius between microrobots can contribute to a significant coverage in a short interval of time. However, due to spatial constraints this is not conserved at medium and long times. From Fig. 3b one can notice that the swarms exhibit a non-monotonic behavior in terms of number of distinct covered/visited spheres as a function of the interaction radius. This has significant implications on the design of interacting micro-robotic swarms. Indeed, by designing the interacting micro-robots independently of the spatial constraints can lead to a fast search in the early times but also to a saturation point in terms of the region visited by the micro-robots and this can lead to inefficient detection of any abnormal behavior. For the current experimental setup, one can safely assume that an interaction
radius of 20 hops offers the best detection of the abnormal behavior. Nevertheless, in practical cases the interaction radius need to be carefully set as a function of the space in which the micro-robots move while considering the characteristics of the movement process (i.e., the directional transition probabilities). Besides coverage, it is also essential to quantify how much time it takes for a group of micro-robots to reach a certain targeted region; this would be a hitting time metric [5]. Consequently, we measured the hitting times for all microrobots to reach their predefined destinations for four interaction radius: 1 hop (1 μm), 5 hops (5 μm), 10 hops (10 μm), 15 hops (15 μm), respectively. Fig. 4 shows the probability of hitting time to exceed a certain threshold as a function the four interaction radius. The hitting time probabilities were also obtained from 1000 independent simulation runs. One can observe that by considering larger interaction radii, the probability of hitting time to exceed a given threshold tends diverge from the Gaussian law and is better approximated by long tail type of distributions such as Fisk, Student’s t and generalized extreme value. V. CONCLUSIONS AND FUTURE WORK In this paper, we have set forth a non-equilibrium statistical physics description of micro-robotic swarms with application to abnormal behavior detection and drug delivery. Unlike previous swarm models, we have developed a mathematical formalism which takes into consideration the volume exclusion rules and micro-robot-to-micro-robot interactions. To get more intuition about the implications of spatial dimensions of microrobots and their interactions, we have presented a comprehensive simulation investigation. Our results show that volume exclusion rules and interaction radii have significant implications on the probability of detecting abnormal behavior (e.g., cancer cells, fungal infections) at both short and long times. In addition, volume exclusion rules and interaction radii also affect the distribution of drug delivery (hitting) times. Consequently, when designing dense networks of microrobotic swarms it is essential to take into account the dimensions of micro-robots, their volume exclusion rule, and the nature of their interactions. ACKNOWLEDGMENT This work was supported by the US National Science Foundation (NSF) under Grant CPS-1135850.
a)
b)
c)
d)
Fig. 4 Probability (Hitting times> threshold) for various interaction distances among micro-robots within a 1cm×1cm×1cm region. a) 1 hop interaction. b) 5 hops interaction. c) 10 hops interaction. d) 15 hops interaction. In scenarios c) and d), the distribution of hitting times is better approximated by means of the generalized extreme value, Fisk and Student’s t-distributions rather than by means of the Gaussian law.
References
R. Axelrod, The Complexity of Cooperation: Agent-Based Models of Competition and Collaboration, Princeton University Press, 1997. [2] R. E. Baker, C. A. Yates and R. Erban, “From Microscopic to Macroscopic Descriptions of Cell Migration on Growing Domains,” Bulletin of Mathematical Biology, Vol. 72, No. 3, pp. 719-762, 2010. [3] B. Behkam and M. Sitti, “Characterization of Bacterial Actuation of Micro-Objects,” IEEE Intl. Conf. on Robotics and Automation, 2009. [4] H. C. Berg, Random Walks in Biology, Princeton Univ. Press, 1993. [5] P. Bogdan and R. Marculescu, “Hitting Time Analysis for Fault-Tolerant Communication at Nanoscale in Future Multiprocessor Platforms,” IEEE Trans. on CAD of Integrated Circuits and Systems, vol.30, no.8, 2011 [6] C. M. Breder, “Equations Descriptive of Fish Schools and Other Animal Aggregations,” Ecology, Vol. 35, No. 3, pp. 361-370, Jul., 1954. [7] D. L. Deangelis and L. J. Gross (Editors), Individual-Based Models and Approaches in Ecology: Populations, Communities and Ecosystems, Chapman and Hall. 1992. [8] A. Galstyan, T. Hogg, and K. Lerman, “Modeling and Mathematical Analysis of Swarms of Microscopic Robots,” IEEE Swarm Intelligence Symposium (SIS-2005), Pasadena, CA, 2005. [9] V. Gazi and K.M. Passino, Swarm Stability and Optimization, Springer Technology & Engineering, 2011. [10] S. Gueron and S. A. Levin, “The Dynamics of Group Formation,” Mathematical Biosciences, Vol. 128, Issues 1–2, pp. 243-264, 1995. [11] E. Gul, B. Atakan, and O.B. Akan, “NanoNS: A Nanoscale Network Simulator Framework for Molecular Communications”, Nano Communication Networks Journal (Elsevier), vol. 1, pp. 138-156, 2010. [12] D. Hanahan and R.A. Weinberg, “The Hallmarks of Cancer,” Cell, Vol. 100, Issue 1, pp. 57-70, January 2000. [1]
[13] C. Kipnis and C. Landim, Scaling Limits of Interacting Particle Systems, Springer - Mathematics, 1999. [14] I. Llatser, I. Pascual, N. Garralda, A. Cabellos-Aparicio, and E. Alarcón, “N3Sim: A Simulation Framework for Diffusion-based Molecular Communication”, IEEE TC on Simulation, no. 8, pp. 3-4, March 2011. [15] R. Mahnke, J. Kaupužs and I. Lubashevsky, Physics of Stochastic Processes: How Randomness Acts in Time, Wiley Mathematics, 2009. [16] A.Mogilner and L.Edelstein-Keshet, “A Non-Local Model for a Swarm,” Journal of Mathematical Biology, vol. 38 , pp. 534-570, 1999. [17] J. D. Murray, Mathematical Biology, Springer - Science, 2003. [18] G. Naldi, Mathematical Modeling of Collective Behavior in SocioEconomic and Life Sciences, Springer-Science, 2010. [19] F. Schweitzer and J. D. Farmer, Brownian Agents and Active Particles: Collective Dynamics in the Natural and Social Sciences, Springer, 2007. [20] M. Sitti, “Miniature Devices: Voyage of the microrobots,” Nature, vol. 458, pp. 1121-1122, 30 April 2009. [21] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, “Novel Type of Phase Transition in a System of Self-Driven Particles,” Phys. Rev. Lett. , vol. 75, pp. 1226–1229, 1995. [22] K. Warburton and J. Lazarus, “Tendency-Distance Models of Social Cohesion in Animal Groups,” Journal of Theoretical Biology, 1991. [23] R. Weiss and T.F. Knight, Jr. “Engineered Communications for Microbial Robotics,” 6th Intl. Workshop on DNA-Based Computers: DNA Computing, A. Condon, G. Rozenberg (Eds.) Springer, 1-16, 2000. [24] C. Xue and H. G. Othmer, “Multiscale Models of Taxis-Driven Patterning in Bacterial Populations,” SIAM Journal on Applied Mathematics, vol. 70, pp. 133-167, 2009. [25] C.A. Yates, R.E. Baker, R. Erban, and P.K. Maini, “Refining SelfPropelled Particle Models for Collective Behaviour,” Canadian Journal of Applied Mathematics, vol. 18, pp. 299-350, 2010.