University of Southern California Dissertations and Theses

Mathematical characterizations of microbial communities: analysis and implications

(USC Thesis Other)

Mathematical characterizations of microbial communities: analysis and implications

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Mathematical Characterizations of Microbial Communities; Analysis and Implications by Hana Koorehdavoudi A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSIPHY (MECHANICAL ENGINEERING) August 2017 Copyright 2017 Hana Koorehdavoudi ii DEDICATION To my family for their everlasting love and support iii Acknowledgments I would like to express my special appreciation and thanks to my thesis advisors and mentors, professor Paul Bogdan and professor Henryk Flashner, for their motivation, support and invaluable guidance through my Ph.D. studies at University of Southern California. I would like to thank them for encouraging my research and allowing me to grow as an independent thinker and research scientist. I would like to thank Professor Paul Bogdan, who inspired me to think outside the box and dream about novel solutions to various research problems. I would like to show my gratitude to Professor Henryk Flashner for all of his priceless advices during the past years. His support during these years made it possible for me to reach this stage in my life. I also would like to extend my sincerest thanks to my defense committee member, professor Geoffrey Spedding for his great comments and suggestions to enhance the quality of my research and his great support. A special thanks to my family. Words cannot express how grateful I am to my parents, Ziba and Rahmat, and my brothers, Kasra and Kian. Their unconditional love, support, encouragement and advice through life helped me to strive towards my goals. Specially, I would like to express my deepest appreciation to my parents, who has been on my side in all of my tough times and sleepless nights through my doctoral studies, while being physically thousands of miles away from me. Without them I could not achieve any of my dreams in life. Last but not least, I have been privileged to benefit from the guidance of professors and researchers from Carnegie Mellon University. I am especially grateful to professor Metin Sitti, professor Radu Marculescu for their help and guidance. iv Contents 1 INTRODUCTION ...............................................................................................................................1 2 BACTERIUM INTERNAL DYNAMICS .........................................................................................5 2.1 INTRODUCTION ...............................................................................................................................5 2.2 CHEMORECEPTOR DYNAMICS PERFORMANCE IN CHEMOTAXIS PATHWAY ..................................6 2.3 TWO DIFFERENT MODELS FOR COOPERATIVITY OF CHEMORECEPTORS ......................................8 2.3.1 Ising Model ............................................................................................................................8 2.3.2 Monod-Wyman-Changeux (MWC) Model ...........................................................................9 2.4 COMPARING TWO DIFFERENT METHYLATION KINETICS MODELS .............................................10 2.4.1 Fine-Tuned Model ...............................................................................................................10 2.4.2 Barkai and Leibler Model ....................................................................................................11 2.5 COMBINATION OF THE MENTIONED MODELS TO SIMULATE STANDARD MODEL OF E. COLI. ...13 2.5.1 Cooperative Chemoreceptor Cluster ....................................................................................13 2.5.2 Phosphorylation Relay and Flagellar Motor ........................................................................15 2.5.3 Bacteria Transport Model: ...................................................................................................15 3 MATHEMATICAL CHARACTERIZATION OF SINGLE BACTERIUM MOTION ............17 3.1 INTRODUCTION .............................................................................................................................17 3.2 SIMULATIONS ...............................................................................................................................18 3.3 EXPERIMENTS ...............................................................................................................................21 3.3.1 Bacteria Culture ...................................................................................................................21 3.3.2 Experimental Setup ..............................................................................................................21 3.3.3 Experimental Conditions .....................................................................................................22 3.3.4 Imaging and Tracking ..........................................................................................................22 3.4 INVESTIGATION OF SIMULATED E. COLI SWIMMING DYNAMICS ................................................23 3.4.1 Introduction on Ergodicity Test ...........................................................................................23 3.4.2 Exploring the Second Order Moment of Bacterium Motion with the Ergodicity Test .......23 3.4.3 Introduction on Nonlinearity Tests ......................................................................................28 3.4.4 Exploring the Third Order Moment of Bacterium Motion with Nonlinearity Tests ...........30 3.4.5 Introduction on Multi-fractal Detrended Fluctuation Analysis (MFDFA) ..........................33 v 3.4.6 Exploring Higher Order Moment of Motion Using MFDFA Analysis ...............................36 3.5 INVESTIGATION OF SIMULATED S. MARCESCENS SWIMMING DYNAMICS ...................................39 3.6 INVESTIGATION OF REAL S. MARCESCENS SWIMMING DYNAMICS THROUGH IN VITRO EXPERIMENTS. ..........................................................................................................................................40 3.7 SUMMARY ....................................................................................................................................41 4 MATHEMATICAL CHARACTERIZATION OF BACTERIA COLLECTIVE MOTION ....43 4.1 INTRODUCTION .............................................................................................................................43 4.2 ENERGY LANDSCAPE ...................................................................................................................44 4.3 ALGORITHM TO CONSTRUCT THE ENERGY LANDSCAPE OF A COLLECTIVE MOTION .................45 4.4 CONSTRUCTING THE ENERGY LANDSCAPE OF A SIMULATED COLLECTIVE MOTION .................46 4.5 INTERPRETATION OF MISSING INFORMATION FOR A COLLECTIVE MOTION ...............................50 4.6 QUANTIFYING MISSING INFORMATION FOR A SIMULATED COLLECTIVE MOTION .....................51 4.7 INTERPRETATION OF EMERGENCE, SELF-ORGANIZATION AND COMPLEXITY FOR A COLLECTIVE MOTION ....................................................................................................................................................52 4.8 COMPLEXITY ANALYSIS FOR THREE NATURAL COLLECTIVE MOTIONS .....................................54 5 MATHEMATICAL ANALYSIS OF BACTERIA SWARMS FOR FIGHTING INFECTIONS: A CASE STUDY ON EBOLA TREATMENT APPROACHES ...........................................................62 5.1 INTRODUCTION .............................................................................................................................62 5.2 EBOLA VIRUS DYNAMICS WHILE INFECTING THE HUMAN CELLS ..............................................63 5.3 MATHEMATICAL MODEL OVERVIEW ..........................................................................................64 5.4 VIRUS INFECTION MODEL RESISTANCE TO CHEMOTHERAPY .....................................................65 5.5 VIRUS INFECTION MODEL RESISTANCE TO GENE THERAPY .......................................................66 5.6 VIRUS INFECTION MODEL RESISTANCE TO BACTERIA THERAPY ...............................................67 5.7 SYSTEM IDENTIFICATION .............................................................................................................68 5.8 SIMULATION RESULTS AND ANALYSIS ........................................................................................70 6 CONCLUSION AND FUTURE RESEARCH DIRECTIONS ......................................................74 6.1 MAJOR CONTRIBUTIONS OF THIS DISSERTATION ........................................................................74 6.2 FUTURE RESEARCH DIRECTION ...................................................................................................75 7 BIBLIOGRAPHY ..............................................................................................................................76 vi List of Tables Table 3.1 Comparing TAMSD and EMSD for different densities of cases 1 for simulated E. coli bacterium trajectory ......................................................................................................................................24 Table 3.2 Slope of line fitted to TAMSD plot of a single bacterium ...........................................................26 Table 3.3 Summary of linearity tests on a single bacterium trajectory from simulated E. coli bacteria. These result are based on different tests on the bacterium trajectory from simulated E. coli with density of 8×10 2 (Bacteria/cm3) in a cubic environment with size 5mm×5mm×5mm without chemoattractant in the environment. Number of variables in all the tests are 3 as the bacterium trajectory is in 3D environment, and the sample size is 1001. All the tests results are in agreement that the bacteria motion does not have normal distribution; hence, bacterium has nonlinear motion. .........................................................................................................................................33 Table 5.1 Parameters and variables of Ebola model ....................................................................................69 vii List of Figures Figure 2.1 Bacterium senses chemoattractant in the environment, sends a signal from receptors through its chemical pathway to the flagellar motor. Based on this, a bacterium is able to control the rotation direction of flagella. This process helps the bacteria to recognize and respond to surrounding chemical signals. ..........................................................................................................................6 Figure 2.2 E. coli chemotaxis signaling pathway ..........................................................................................7 Figure 2.3 Models of cooperativity of receptors. a) All-or-none MWC model. In this model, the chemoreceptors group into smaller sub-clusters. All the receptors in each cluster are coupled and have the same state. b) Two-state Ising model. In this model, the cooper ..........................10 Figure 2.4 Two simplified models for methylation and demethylation process in E. coli bacteria, which capture exact adaptation property. a) Fine-Tuned model. In this model, the methylated receptors are active; while the unmethylated receptors are inactive. Therefore, exact adaptation depends on proportion of different biochemical parameters. b) Barkai-Leibler model. In this model, unmethylated receptors are in inactive state; while the methylated receptors can transit between active and inactive states. Therefore, exact adaptation is robust and it is valid for a wide range of biochemical parameters .........................................................................................................13 Figure 2.5 Simulation results for chemoreceptor activity in E. coli showing precise adaptation. The figure is obtained by averaging simulation results over 370 clusters of receptors [60]. ......................14 Figure 3.1 Simulation configuration of E. coli in BNSim (see Methods) and single particle tracking of S. marcescens in vitro experiments. (a) In case 1 for simulation of E. coli in BNSim, a cubic environment with dimensions 5000x5000x5000 µm 3 has been considered. The whole environment subdivided into 10 6 smaller cubes each with size of 50x50x50 µm 3 . We consider a linear gradient (10 -4 mM/µm) of chemoattractant (L-aspartate) only in the y-direction and no gradient in other directions. (b) In case 2 for simulation of E. coli in BNSim, we consider the same configuration as case 1 for bacteria and environment. We only changed the gradient of chemoattractant in the environment. In this case we consider a linear gradient of chemoattractant (L-aspartate) from the injection location to the targets, meaning there is linear gradient in 3 different dimensions. (c) Single particle tracking of S. marcescens in an in vitro environment viii with dimensions of 10000x500x150 µm 3 with linear gradient (10 -4 mM/µm) of chemoattractant (L-aspartate) in the second direction (y-direction). ...................................................................20 Figure 3.2 Mean square displacement (MSD) plots for simulated E. coli in BNSim. (a) Time averaged MSD of one E. coli bacterium motion simulated in BNSim as a function of lag time Δ (δ2(Δ,T) = 1T-Δ0T-Δr (t+Δ) -r (t)2dt) for a density of 8x10 4 (bacteria/cm 3 ) without chemoattractant and obstacles. This plot shows superdiffusion behavior (MSD (t) = kαtα,1 < α < 2). (b) Ensemble averaged MSD of one E. coli motion simulated in BNSim vs. time t calculated from the following formula (MSD (t) = < |r (t) -r (0)|2 >) over 100 trajectories for a density of 8x10 4 (bacteria/cm 3 ) without chemoattractant and obstacles. This plot also shows superdiffusion. The observed differences between these two plots (a and b) show an ergodicity breaking in the system and demonstrate the nonergodicity of bacteria dynamics. (c) The TAMSD of the center of whole swarm for E. coli motion simulated in BNSim shows the effect of population density and chemoattractant on center of swarm motion (case 2). (d) The TAMSD of one E. coli trajectory simulated in BNSim for different cases show the effect of chemoattractant and density on a single bacterium motion (case 2). By increasing the population density, we observe a phase transient in bacteria behavior from superdiffusion to normal diffusion and then subdiffusion (case 2). ...................................................................................25 Figure 3.3 Effects of time, bacteria population density and chemoattractant on E. coli motion simulated in BNSim (a) Effect of time on bacteria motion. Log-log plot of P(r,τ) as a function of distance r for different lag times τ = 100 s,500 s,1000 s. This result is for 8x10 2 bacteria/cm 3 in an environment without chemoattractant and obstacles. This plot shows that the bacteria move more with time and get further away from their initial locations. The peak of the plot decreases more (case 1). (b) Effect of density on bacteria motion. Log-log plot of P(r,τ) for lag time τ = 1000 s compared for different bacteria (E. coli) density simulated in BNSim in an environment with chemoattractant gradients and without obstacles. By increasing the bacteria population, their motion will be restricted and they are able to move less freely in the environment (case 1). (c) Effect of chemoattractant on the bacteria motion. Log-log plot of P(r,τ) for lag time τ = 1000 s compared for the case with and without chemoattractant in the environment for the population of 8x10 2 bacteria/cm 3 . By adding chemoattractant in the environment the bacteria tend to move more and oscillate less so they get further away from their initial conditions in the direction of increasing gradient of chemoattractant (case 1). ....................................................27 Figure 3.4 Skewness plots for some of the simulation E. coli trajectories. a) 8x10 2 (bacteria/cm 3 ) in the cubic environment, with/out linear gradient of chemoattractant in y direction, with/out obstacles on the bacteria way. b) 8x10 3 (bacteria/cm 3 ) in the cubic environment, with/out linear gradient ix of chemoattractant in y direction, with/out obstacles on the bacteria way. c) 8x10 4 (bacteria/cm 3 ) in the cubic environment, with/out linear gradient of chemoattractant in y direction, with/out obstacles on the bacteria way. In all the cases skewness is unequal to zero meaning bacterium displacement does not have normal distribution; therefore, bacterium motion is nonlinear. ....32 Figure 3.5 Multi-fractal analysis of motion for both simulated E. coli in BNSim and S. marcescens in vitro experiments (a) Multi-fractal spectrum of the simulated E. coli position in the direction of linear gradient of chemoattractant for different cases with and without chemoattractant and different bacteria densities (case 2). (b) Generalized Hurst exponent of the simulated E. coli position in the direction of linear gradient of chemoattractant for different cases with and without chemoattractant and different bacteria density (case 2). (c) Multi-fractal spectrum of the simulated E. coli position in direction of linear gradient of chemoattractant with a density of 8x10 4 bacteria/cm 3 for different cases with and without chemoattractant and obstacles in the environment (case 1). (d) Generalized Hurst exponent of the simulated E. coli position in direction of linear gradient of chemoattractant with a density of 8x10 4 bacteria/cm 3 for different cases with and without chemoattractant and obstacles in the environment (case 1). (e) Multi- fractal spectrum of the experiments in vitro on S. marcescens position in the direction of linear gradient of chemoattractant for different cases with and without chemoattractant and different bacteria densities. (f) Generalized Hurst exponent of the experiments in vitro on S. marcescens position in the direction of linear gradient of chemoattractant for different cases with and without chemoattractant and different bacteria densities. (g) Multi-fractal spectrum of simulated S. marcescens position in the direction of linear gradient of chemoattractant for different cases with and without chemoattractant and different bacteria densities. (h) Generalized Hurst exponent of simulated S. marcescens position in the direction of linear gradient of chemoattractant for different cases with and without chemoattractant and different bacteria densities. .....................................................................................................................................38 Figure 3.6 Skewness plots for simulated s. marcescenns trajectories. Skewness plots for different bacteria density cases in simulation, with/out linear gradient of chemoattractants in y direction. In all the cases skewness is unequal to zero meaning bacterium displacement does not have normal distribution; therefore, bacterium motion is nonlinear. .............................................................39 Figure 3.7 Skewness plots for experimental S. marcescens trajectories. Skewness plots for different bacteria density cases in vivo environment, with/out linear gradient of chemoattractant in y direction. In all the cases skewness is unequal to zero meaning bacterium displacement does not have normal distribution; therefore, bacterium motion is nonlinear ..............................................................40 Figure 4.1 Energy landscape of a complex system with one degree of freedom. ........................................45 x Figure 4.2 Schematic description of the main steps for building the energy landscape for a group of N agents moving in a three-dimensional space. (a) First, we subdivide the trajectories of all agents in the group to equal sub-intervals centered at time tc with a time window of [tc-∆2,tc+∆2], where ∆ is the predefined time scale. Next, we estimate the three-dimensional probability distribution function of the motion of the group for each sub-interval. (b) We use the Kantorovich metric to cluster these sub-interval time series based on their similarities in the probability distribution function. Each cluster of sub-intervals can be interpreted as a state for the collective motion. (c) In the last step, we estimate the transition probability matrix among the identified states of the collective motion. .......................................................................................................................46 Figure 4.3 Different zones of interaction around each individual in a group of agents moving in three- dimensional space in a model proposed by Couzin and his coworkers [146]: Zone of repulsion, zone of orientation and zone of attraction. There is a blind zone behind each individual, in which they do not sense and react towards other neighbors in that area. .............................................48 Figure 4.4 Various collective patterns of a simulated model of a group of agents moving in a three- dimensional space. (a) Torus: Individuals rotate around a center point within an empty space. (b) Swarm: Individuals show attraction and repulsion behavior between themselves and there is no orientation behavior consequently no parallel motion. (c) Dynamic parallel group: Individuals align with each other and make the group more motile compared to two previous cases. (d) Highly parallel group: Individuals are in a highly aligned arrangement and the group is more motile compared to previous cases. ..............................................................................49 Figure 4.5 Transition probabilities among the states identified in different collective behaviors of the simulated agent-based model. (a) Torus, the plot shows the transition probability between different states in this collective behavior. (b) Swarm, the group of agents has the highest number of states in this collective behavior and the landscape has more spikes and is less smooth compared to the other cases. (c) Dynamic parallel group, the transition probability looks similar to the torus phase, this similarity is due to preference of individuals to align their motion parallel to their neighbors. (d) Highly parallel group, the group has the lowest number of possible states in this phase and the landscape is less spiky due to high preference of individuals to move parallel with respect to each other. ............................................................................................50 Figure 4.6 Quantifying the missing information of the entire simulated agent-based model for various interaction rules. (a) We quantify the missing information from the dynamics of a group of agents considering different interaction rules which causes various collective behavior in the group while considering the same initial condition for the agents position. This plot shows the transition from swarm phase to torus, dynamic parallel group and then highly parallel group xi and the fact that the missing information is decreasing due to an increase in the internal order of the group. (b) The quantified missing information extracted from the group dynamics when the population consists of 100 individuals for varying value of the radius of zone of orientation while radii of zone of repulsion and attraction is fixed (red line) and in other case varying value of the radius of zone of attraction while radii of zone of repulsion and orientation is fixed (blue line). ...........................................................................................................................................52 Figure 4.7 Transition probability matrix and complexity analysis for a collective group of bacteria. (a) Transition probabilities among the possible states for a collective group of 9 bacteria selected from a population density of 10 8 bacteria/cm 3 moving in an environment without chemoattractant gradient. (b) Complexity analysis for different states compared to the first identified state with lowest energy level for a group of 9 bacteria selected from a population density of 10 8 bacteria/cm 3 moving in an environment without chemoattractant gradient. This plot shows the level of change in missing information when the collective motion leaves each identified state to evolve to a new state. It also demonstrates the relative emergence and relative self-organization and relative complexity of the collective motion when evolving from any of the identified state to the first and most probable state. (c) Transition probabilities between the possible states for a group of 9 bacteria selected from a population density of 10 8 bacteria/cm 3 moving in an environment with chemoattractant gradient. (d) Complexity analysis for different states compared to the first identified state with lowest energy level for a group of 9 bacteria selected from a population density of 10 8 bacteria/cm 3 moving in an environment with chemoattractant gradient. ...........................................................................................................55 Figure 4.8 Level of change in missing information through the evolution of the dynamic of collective group from state 1 to other states. ........................................................................................................56 Figure 4.9 Transition probability matrix and complexity analysis for a collective group of pigeons. (a) Transition probabilities between the possible states for a group of 9 pigeons in free flight. (b) Complexity analysis for different states compared to the first identified state with lowest energy level for a group of 9 pigeons in free flight. (c) Transition probabilities between the possible states for a group of 8 pigeons in home flight. (d) Complexity analysis for different states compared to the first identified state with lowest energy level for a group of 8 pigeons in home flight. ..........................................................................................................................................57 Figure 4.10 Transition probability matrix and complexity analysis for a collective group of ant. (a) Transition probabilities between the possible states for a group of ants. (b) Complexity analysis for different states compared to the first identified state with lowest energy level for a group of ants. ............................................................................................................................................58 xii Figure 4.11 The proposed algorithm can be integrated into an engineered framework which could help to set the parameters representing the dynamics of the agents and their interactions to achieve a certain degree of emergence, self-organization and complexity. ..............................................61 Figure 5.1 Diagram of Ebola infection. Cell level interaction of Ebola virus with the human immune system T cells. Virus inters human T cells and infects them. Bursting infected cells release copies of virus............................................................................................................................................64 Figure 5.2 Diagram of gene therapy for Ebola infection. Cell level dynamics interactions of Ebola virus and the human immune system T cells during gene therapy. Ebola virus can enter both T cells and P cells and infect them. ..............................................................................................................66 Figure 5.3 Fluctuation of different species in the cell level network of Ebola infection. a) Population dynamics of healthy and infected T cells during Ebola infection. Healthy T cells population reaches less than 10% of its initial population after seven days. b) Growth of virus in the first ten days of infection. ..................................................................................................................70 Figure 5.4 Population fluctuation of cell level network of immune system T cells and Ebola virus during chemotherapy. Symbol ∆ represents the treatment period and P is the efficiency of the drug during treatment. For different cases of drug efficiency with various treatment period the T cell population reaches 20% or less than its initial population in the first 30 days of Ebola infection. ....................................................................................................................................................71 Figure 5.5 Fluctuation of different species in the cell level network of Ebola infection during gene therapy. a) Population dynamics of healthy and infected T cells and also healthy protected P cells during Ebola infection. b) Fluctuation of infected protected P cells c) Growth of virus under gene therapy in the first forty days of infection. ................................................................................71 Figure 5.6 Population fluctuations of different species in the network of immune system T cells, Protected P cells and Ebola virus during gene therapy. T 0 is the initial population of unprotected CD4+T cells, and P 0 is the initial population of protected cells. T cells include both protected and unprotected cells population. Based on these results, decreasing the ratio of initial unprotected to protected T cells results in less reduction in T cell population, which shows higher efficiency of gene therapy. ..........................................................................................................................72 Figure 5.7 Fluctuation of different species in the cell level network of Ebola infection during bacteria therapy. a) Population dynamics of healthy cells during Ebola infection. b) Fluctuation of infected T cells c) Growth of virus under bacteria therapy in the first forty days of infection. 73 Figure 5.8 Comparison of three different therapeutic strategies. Bacteria therapy has less damage on uninfected cell population and it will cure the infection in the shortest amount of time. ..........73 1 1 Introduction The future potential of using bacteria for therapeutic purposes and regenerative medicine makes the dynamics of such micro-swimmers highly attractive to study [1]. As an example, Felfoul and his coworkers show that magnetotic bacteria can be used to transport drug-loaded nanoliposomes into oxygen-depleted hypoxic region of the tumor [2]. Considering such applications, significant research efforts have been focused on analyzing and modeling bacteria swimming dynamics. Broadly speaking, the mathematical models used to describe the bacterial swimming dynamics can be classified into two categories. The first category is based on a microscopic (i.e. cell-level) view of the bacterial swimming. In this category, a set of equations describes the state of a single agent [3]–[7]. The second category provides a macroscopic (i.e. population-level) view via continuum-based partial differential equations that capture the dynamics of population density over space and time, without considering the intracellular characteristics directly [8]– [14]. As a brief overview, models presented by Cates [15] highlight that bacterial dynamics does not always obey detailed balance which means it is a biased diffusion process depending on the environmental conditions. Along this direction, Schnitzer [16] uses the Smoluchowski equation to describe the biased random walk of the bacterium Escherichia coli (E. coli) during chemotaxis to search for food. Croze and his coworkers [17] study experimentally and theoretically the effect of concentration of soft agar on chemotaxis of E. coli as an example of environmental condition affecting this biased random walk. Gennes [18] derives the average run length traveled by E. coli during one counter-clockwise interval to present a detailed description of the motion taking place during one run interval. Ariel and his coworkers show that the bacteria preform super-diffusion during swarming on a surface [19]. Indeed, despite the research and progress in this area, the researchers in this field are still trying to define more appropriate mathematical tools for describing bacteria motion and predicting their behavior both in vivo and in vitro, as the models proposed to date still remain mostly imprecise. This is because these models are developed under simple assumptions (e.g., bacteria swim in a free space, chemoattractant diffuse uniformly, etc.) and do not take into account realistic conditions like volume exclusion (i.e., no two bacteria can occupy the same space at the same time), chemical interactions among bacteria (autochemotaxis 1 ), 1 Autochemotaxis happens when bacteria excrete the converted substrate succinate molecules in the environment into chemoattractant aspartate molecules. Once the succinate gets depleted, bacteria start consuming the chemoattractant 2 obstacles, and heterogeneous distribution of chemoattractant in the environment. Consequently, new mathematical models are needed to take into account the effects of such realistic conditions upon individual bacterium in a population. Developing such microscopic mathematical characterizations can also lead to modeling the bacteria population dynamics more accurately at the macroscopic scale. Furthermore, developing such models can help with a more accurate prediction and control of the bacterial swarm dynamics inside the body to match the patient’s physiological conditions. Toward this end, in this thesis, we identify the main characteristics of a single bacterium motion, which could have a fundamental impact on mathematical modeling of both single bacterium motion and swarm dynamics. We analyze the motion trajectories of E. coli and Serratia marcescens (S. marcescens); these are two types of bacteria which have very similar swimming dynamics, behavior, and chemotaxis [20], [21]. Extracting single bacterium trajectories of motion from computer simulations for E. coli and S. marcescens and applying single-cell visual tracking methods to real experiments on S. marcescens reveal the main mathematical characteristics of these bacteria motion. To elucidate the complexity of the bacterial dynamics in cell-level, we investigate the individual trajectories of simulated E. coli and S. marcescens beside trajectories of real S. marcescens at multiple scales in space and time. This way, we show for the first time that their motion exhibits a nonlinear and multi-fractal behavior. The observed nonlinearity in the individual based bacterium motion can be attributed to the nonlinear dependencies within the bacterium chemical pathway known as internal dynamics. Moreover, the multi-fractal feature means that bacterial motion is self-similar, which is a concept dating back to the pioneering work of Kolmogorov to explain the chaotic behavior [22], [23]. The multi-fractal formalism offers mathematical techniques to capture this nonlinear self-similarity in their motion. Overall, our new results suggest that multi-fractal modeling can be a new modeling paradigm for bacterial swarms in microscopic level. In the next step, we study the collective motion of bacteria in macroscopic level and quantify the complexity of their group dynamics. Their collective motion is an example of a complex system. Generally speaking, a complex system refers to a system in which there is a lack of precise relation between the system outcomes and the original causes of those outcomes, which causes the system to have unpredictable and nonlinear dynamics. The complexity in such a system is due to the tight coupling between the components of the system, which makes it impossible to analyze the components individually and isolated from the rest of the system[24]. This close coupling and interactions between the units of the complex system cause recognizable collective behavior at a larger scale [25]. aspartate excreted by themselves and they may get some information about the nutrition conditions and population of bacteria around them through this process [109], [124]–[126]. 3 A group of agents or animals moving collectively is an example of a complex system. The collective characteristics between the agents of a group are regulated by behavioral tendencies, as well as short-range and long-range interactions among them. In such a complex system, the group with identical agents’ behavior evolves through different states (i.e., spatio-temporal arrangement/configuration of the agents moving in a collective group formation) [26], [27]. We can encode the dynamics (evolution) of the group among different states by constructing a free energy landscape representation. Different factors like the number of members, internal capabilities of the individuals (e.g., sensitivity to neighbors, motion speed of individuals, computational/processing capabilities of agents) and external properties (e.g., environmental and boundary condition) influence the overall collective behavior and the free energy landscape; further, these factors can contribute to various phase transitions in the group structure among several possible states in the corresponding free energy landscape. Despite significant research and progress in studying natural [28]–[36] and engineered [37]–[49] collective motion, the scientists in the field still try to quantify the collective dynamical state in a swarm motion and predict the transition between these dynamically stable states. Toward this end, we developed a new approach, which extracts the dominant states of spatial formation and structures for a swarm. Our framework presents the free energy landscape [50]–[62] of the states of the swarm and also the transitions between these states. In this approach, we are able to distinguish between absorbing and transition states. Absorbing states are the ones in which the swarm prefers to stay for a relatively longer time. In contrast, the swarm forms into a transition states while transitioning from one absorbing state to another absorbing state. Moreover, based on the energy landscape, we are able to quantify the missing information related to each of the swarm states and also the emergence, self-organization and complexity of the structural formation of each state of swarm with respect to other states. The measures we adopt are based on Shannon’s definition of information [63]. Our approach enables a mathematical quantification of biological collective motion complexity, which allows us to recognize and differentiate among various possible states based on their relative energy level and complexity measures. Our new results both in microscopic and macroscopic level can help us better investigate new bacteria- based approaches for a variety of purposes like diagnostic and targeted drug delivery. This new results can help to improve the prediction, optimization and control of bacterial population dynamics in crowded environments. This can prove beneficial for constructing mathematical models of the swarm of bacteria propelled micro-robots and quantify their performance. The remaining of this thesis is organized as follows: In the second chapter, we briefly overview a bacterium internal dynamics including chemoreceptors, chemical pathway and flagella motor. We present different models for chemoreceptors that capture the sensitivity and adaptively of bacterium to 4 environmental condition. In the third chapter, we present our analysis which reveals the main characteristics of single bacterium motion. In this chapter, we show the bacterium motion is nonergodic, nonlinear and multi-fractal. In chapter four, we present our framework to extract the possible states in swarm formation and the strategy to build the corresponding energy landscape for transitions between these possible states. In this chapter, we also define missing information for the possible states of a swarm structure and measure emergence, self-organization, and quantify the complexity of a swarm based on those. In chapter five, we present a case study for using bacteria for therapeutic purposes. More precisely, we develop a mathematical model for the interactions of Ebola virus with human immune system. Based on this mathematical model, we propose bacteria therapy as a potential approach to cure infection like Ebola based on our mathematical model. 5 2 Bacterium Internal Dynamics 2.1 Introduction Bacteria can sense particular chemicals in their environment and adapt their behavior and motion with respect to the gradient of the chemical concentration. This phenomenon is referred as bacterial chemotaxis. The chemicals that affect the bacteria motion are known as ligands, which can be categorized into two groups: attractants and repellents, respectively. Bacteria tend to move toward attractants, while move away from repellents. Attractants and repellent molecules are recognized by special protein called chemoreceptor complex. Chemoreceptor is a protein molecule in the surface of a cell, which receives the chemical signals from outside the cell and based on that it influences the internal activity of the cell. In other words, chemoreceptors pass information from the outside to the inside of the cell and help the cell to recognize and respond to surrounding chemical signals. The signals received from receptors are translated by signaling chemical pathway and affect the performance of bacteria motor. Motor controls the frequency of changes in rotation direction of flagella. When all the flagella rotate counterclockwise (CCW), the bacterium runs. When one or some of the flagella change their rotation direction to clockwise (CW), the bacterium stops running and keep tumbling in the same location to choose a new direction for its next running phase. Considering bacteria as our system, the input to the system is the concentration of attractant in the environment, and output is the probability of the motor turns CW. This probability regulates the cell’s tumbling frequency. Figure 2.1 exhibits major players in bacterium chemotaxis. Two essential characteristics of bacteria chemotaxis pathway are robust adaptation and signal amplification, which results in ultra-sensitivity of bacteria to changes in attractant concentration [64]–[67]. These two main characteristics are due to receptor dynamics and performance, which are our main focus in this chapter. In this chapter we review the performance of chemical pathway inside bacteria and present the related models to simulate standard model of E. coli. The remaining of this chapter is organized as follows: In section 2.2, we briefly present the chemoreceptor dynamics and chemical pathway performance. In section 2.3, we explain two different models for cooperativity of chemoreceptors which captures signal amplification characteristic in bacteria: Ising model and Monod-Wyman-Changeux model. Next, we compare two different models for methylation kinetics: fine-tuned model and Barkai and Leibler model. We explain in detail which model performs robust adaptation to the changes in environmental condition. 6 In the last section of this chapter, we simulate standard model of E. coli using a combination of the mentioned models. Figure 2.1 Bacterium senses chemoattractant in the environment, sends a signal from receptors through its chemical pathway to the flagellar motor. Based on this, a bacterium is able to control the rotation direction of flagella. This process helps the bacteria to recognize and respond to surrounding chemical signals. 2.2 Chemoreceptor Dynamics Performance in Chemotaxis Pathway Bacterial chemotaxis signaling pathway consists of transmembrain methyl-accepting chemotaxis protein (MCP) receptors and six different proteins known as CheA, CheB, CheR, CheW, CheY, and CheZ (Figure 2.2). Bacteria identify the chemical signals in their environment by trans-membrane chemoreceptors. Chemoreceptors consist of MCPS and an MCP-like protein. MCPs are linked to cytoplasmic Histidine Protein Kinase (HPK) named CheA through an adaptor protein named CheW. The chemoreceptors form homodimers. This type of homodimer is a macromolecular complex formed by two identical MCPs. Each homodimer can bind with one molecule of ligand. We refer chemoreceptor to a homodimer of MCP [66], [68]. When a ligand binds to a receptor, it changes the activity state of the receptor and protein CheA. The activated CheA gains a phosphate group (PO 4 ) from autophosphorylation and transfer it to either one of two Response Regulator (RR) named CheB and CheY and create CheB-P and CheY-P respectively. The process of adding phosphate group to a protein is called phosphorylation. Phosphorylation is a way for cell to transfer information among signaling proteins. When CheY phosphorylates, CheY-P binds the switch protein FliM on the flagellar motor. Consequently, motor changes its rotate from CCW to CW direction, which causes the bacterium to go from run to tumble mode. Increasing the concentration of CheY-P results in higher tumbling frequency meaning that the bacteria continue swimming towards attractant source [69]. Protein CheZ helps to dephosphorylate the CheY-P and returns the bacterium to run mode [70]. 7 When attractant binds to the receptor, it decreases the autophosphorylation activity of CheA. As a result, it decreases the rate of CheY phosphorylation, and CheY-P concentration drops. Hence, the motor rotates in CW direction with a lower probability. This means, binding of attractant to receptor result in less tumbling frequency. On the other hand, repellents increase the activity of CheA, which has a reverse effect on bacteria behavior. In this condition, bacteria have higher tumbling frequency, which helps them to swim away from the repellent [67], [71], [72]. Bacterium has another mechanism in its chemotaxis pathway, which regulates bacteria adaptation process. When attractant binds to receptor, it will decrease receptor activity. In this condition some methyl groups (CH 3 ) will attach to the receptor result in increasing the activity of receptor. This process is known as receptor methylation. The receptor can have between zero to five methyl groups. Increasing the number of methyl groups cause increasing the activity of the receptor. CheB-P and CheR are the enzymes to demethylate and methylate MCPs and control their functions [73], [74]. The methylation and demethylation circuit continuously works and it is independent of bacterium sensing any ligands. This causes the bacterium to adapt to the changes in ligand’s concentration. In section 2.4, different models of methylation dynamics have been presented which explains the adaptation process in more details. Figure 2.2 E. coli chemotaxis signaling pathway 8 2.3 Two Different Models for Cooperativity of Chemoreceptors In bacteria like E. coli, there are thousands of chemoreceptors. The number of chemoreceptors depends on the growth phase and physiological condition of the bacterium [75]. Chemoreceptors form clusters and perform collectively with each other inside the bacterium [76]. This cooperativity between chemoreceptors causes signal amplification in chemotaxis pathway and increase sensitivity of the bacterium [77]–[81]. Two well-known models of chemoreceptor clusters, which capture cooperativity and signal amplification, are two-state Ising model and all-or-none Monod-Wyman- Changeux (MWC) model. 2.3.1 Ising Model Duke and Bray developed two-state Ising model for receptor based on Ising model for ferromagnetism in physics [82], [83]. In this model, Kinase activity (i.e., phosphorylation process 2 ) of a chemoreceptor has two states: active and inactive, respectively [64]. The active state of the receptor is equivalent to up-spin and the inactive state is equivalent to down-spin. The cooperative interactions between nearest neighbor receptors in the cluster can be modeled based on the Ising ferromagnetic spin-spin interaction to have the same conformation for neighbor receptors (Figure 2.3b) [84]–[89]. The state of a receptor in this model can be distinguished by the combination of two variables (a, l). When a = 0, it represents an inactive receptor. In contrast, when a = 1, it represents an active receptor. When l = 0, it denotes the receptor is empty. In contrast, when l = 1, it shows a binding of the receptor with ligand. The P(a,l) represents the probability that the receptor is in any of the following possible states: P(0,0), P(1,0), P(0,1), and P(1,1). In steady state condition, the relations between the probabilities of the four possible states are as follow: 𝑃 0,1 𝑃 0,0 = 𝐿 𝐾 C (2.1) 𝑃(1,1) 𝑃(1,0) = [𝐿] 𝐾 F (2.2) 𝑃(0,0) 𝑃(0,1) = 𝑒 HI J (K) (2.3) 𝑃(𝑎,𝑙) F,N = 1 (2.4) Where 𝐾 F is the ligand dissociation constant for the active receptor and 𝐾 C is the ligand dissociation constant for the inactive receptor. [L] is the ligand concentration. When there is no ligand bounded to the receptor, 𝑓 K (𝑚) is the free energy difference between active and inactive state which is a function of methylation 2 Kinase activity is a process in which an enzyme adds phosphate group chemically to other protein and modifies them. 9 level m. In the presence of ligand, the free energy difference between active and inactive states is presented by 𝑓 Q ([𝐿]). Hence, the total free energy difference between active and inactive states is presented by equation (2.5): Δ𝑓 = 𝑓 K 𝑚 +𝑓 Q ([𝐿]) (2.5) 𝑓 Q 𝐿 = 𝑙𝑛 1+ [𝐿] 𝐾 F 1+ [𝐿] 𝐾 C (2.6) The averaged activity < 𝑎 > can be determined by equation (2.7) < 𝑎 >= 1 1+𝑒 HSI = 𝑒 I J K 1+ 𝐿 𝐾 F 1+ [𝐿] 𝐾 C +𝑒 I J K 1+ 𝐿 𝐾 F (2.7) In this model, degree of cooperativity can be represented by a correlation length which increases with the receptor interaction strength [66]. 2.3.2 Monod-Wyman-Changeux (MWC) Model In all-or-none MWC model, the chemoreceptors group into smaller sub-clusters [90]. All the receptors in each cluster are coupled and have the same state (Figure 2.3a). The receptors in different sub-clusters are isolated and do not have any cooperativity. In this model, N is the size of receptors in sub-clusters, which shows the degree of cooperativity. Hence, the free energy difference between all the active and inactive states in the same cluster is 𝑁×Δ𝑓 = 𝑁×(𝑓 K 𝑚 +𝑓 Q 𝐿 ). Average activity in this model can be determined by equation (2.8) [66], [90]–[92]. < 𝑎 > = 1 1+𝑒 HU×SI = 𝑒 HU×I J K ×(1+ 𝐿 𝐾 C ) U 𝑒 HU×I J K ×(1+ 𝐿 𝐾 C ) U +(1+ 𝐿 𝐾 F ) U (2.8) 10 Figure 2.3 Models of cooperativity of receptors. a) All-or-none MWC model. In this model, the chemoreceptors group into smaller sub-clusters. All the receptors in each cluster are coupled and have the same state. b) Two-state Ising model. In this model, the cooper 2.4 Comparing Two Different Methylation Kinetics Models Experiments show evidence that when the concentration of chemoattractant (e.g., aspartate) suddenly increases, tumbling frequency of the bacterium drops and then gets back to its initial value [93]. This is a sign of adaptation in bacteria. In general, when a sudden change in the environment of bacterial population happens, the bacterium is able to adapt to the change [94]. In E. coli bacteria, Enzymes CheR and CheB catalyze receptor methylation and demethylation process, which is the source of adaptation. We explain two simplified models for methylation and demethylation process in E. coli bacteria, which capture exact adaptation property. In the first model, exact adaptation is fine-tuned meaning it depends on proportion of different biochemical parameters. In contrast, exact adaptation is robust in the second model and it is valid for a wide range of biochemical parameters [67], [95]. 2.4.1 Fine-Tuned Model This receptor methylation dynamics model was proposed by Knox and colleagues [96]. E. coli receptors get methylated by enzyme CheR, and demethylated by enzyme CheB. In this model, the methylated receptors are active, while the unmethylated receptors are inactive (Figure 2.4a). Variable X represents the number of demethylated receptors and X m is the total number of methylated receptors with activity a 0 per methylated receptor. The methylation enzyme CheR works at saturation with rate V R . On the other hand, demethylation enzyme CheB works with Michaelis-Menten kinetics [97] with rate V B and Michaelis- Menten constant K. The rate of change of receptor methylation X m in this model is as follow: 11 𝑑𝑋 K 𝑑𝑡 = 𝑉 Z 𝑅− 𝑉 ] 𝐵𝑋 K 𝐾+𝑋 K (2.9) The concentration of enzymes CheR and CheB is presented by parameters R and B correspondingly. When the E. coli adapts to the new changes in the chemoattractant concentration in the environment, the methylation dynamics reaches a steady state level with dX m /dt = 0. The level of methylated receptor at steady state condition is: 𝑋 K = 𝐾𝑉 Z 𝑅 𝑉 ] 𝐵−𝑉 Z 𝑅 (2.10) Hence, the total steady-state activity of methylated receptors in E. coli is equal to 𝐴 ` = 𝑎 ` 𝑋 K (2.11) The total activity of the receptor A 0 , due to the receptor methylation, regulates the rate of CheY phosphorylation. Receptor methylation increases the rate of CheY phosphorylation, which results in an increase of tumbling rate. Therefore, the total steady state activity A 0 regulates the steady-state tumbling frequency. When there is a sudden increase in the chemoattractant concentration in the environment, the attractant causes receptors to switch to their inactive conformation. Therefore, a 0 drops to a 1 when a 1 <<a 0 . As a result, the total activity drops to a lower level 𝐴 a = 𝑎 a 𝑋 K , 𝐴 a ≪ 𝐴 ` immediately. This causes the tumbling frequency to decrease initially. Meanwhile the methylation-demethylation cycle continues working. Less number of active receptors cause the rate of action of CheB drops from 𝑉 ] to 𝑉 ] ′; accordingly the demethylation rate decreases. Meanwhile, methylation rates work at saturation. Therefore, receptor methylation 𝑋 K begins to increase. After enough time, the receptor methylation reaches a new steady state level 𝑋 K ′. 𝑋 K d = 𝐾𝑉 Z 𝑅 𝑉 ] ′𝐵−𝑉 Z 𝑅 (2.12) New steady state total activity of the receptors 𝐴 e = 𝑎 a 𝑋 K ′. For exact adaptation the following condition should hold 𝐴 e = 𝐴 ` → 𝑎 ` 𝐾𝑉 Z 𝑅 𝑉 ] 𝐵−𝑉 Z 𝑅 = 𝑎 a 𝐾𝑉 Z 𝑅 𝑉 ] ′𝐵−𝑉 Z 𝑅 (2.13) For this equality to hold, 𝑎 a and 𝑉 ] ′ should balance each other. This means the exact adaptation in this model is fine-tuned and only holds for a precise balance of changes in activity and receptor methylation parameters [67]. 2.4.2 Barkai and Leibler Model We describe a simplified version of the model suggested by Naama Barkai and Stanislas Leibler. In contrast to fine-tuned model, Barkai-Leibler model has a robust exact adaptation, meaning this model reaches exact 12 adaptation for a wide range of parameters. In this model, variable X presents the total number of unmethylated receptors, which are in inactive state. The difference of this model with the previous one is that, in this model, the methylated receptors can transit between active and inactive states (Figure 2.4b). Variables 𝑋 K C and 𝑋 K F represent total number of inactive and active methylated receptors correspondingly. The total activity is equal to the total number of active methylated receptors: 𝐴 ` = 𝑋 K F (2.14) In this model, Barkai and Leibler assumed methylation enzyme CheR works at saturation with rate V R , and demethylation enzyme CheB works with Michaelis-Menten kinetics with rate V B and Michaelis-Menten constant K. Enzyme CheB only demethylate the active receptor 𝑋 K F and do not show effect on the inactive methylated receptors 𝑋 K C . Based on this, the rate of change of methylated receptors 𝑋 K C +𝑋 K F is as follow: 𝑑(𝑋 K C +𝑋 K F ) 𝑑𝑡 = 𝑉 Z 𝑅− 𝑉 ] 𝐵𝑋 K F 𝐾+𝑋 K F (2.15) The concentration of enzymes CheR and CheB presented by parameters R and B correspondingly. When the E. coli adapts to the new level of chemoattractant concentration in the environment, the methylation dynamics reaches a steady state level with 𝑑 𝑋 K C +𝑋 K F /𝑑𝑡 = 0. The level of active methylated receptor 𝑋 K F at steady state condition reaches a point where demethylation exactly balances the constant flux of methylation. The steady state activity: 𝐴 ` = 𝑋 K F = 𝐾𝑉 Z R 𝑉 ] 𝐵−𝑉 Z 𝑅 (2.16) Now imagine that there is a sudden increase in the chemoattractant concentration in the environment. This reduces the probability that the methylated receptor is in active state; hence the number of active receptor drops suddenly which results in lower tumbling frequency. Meanwhile, enzyme CheR works at saturation and methylates receptors with the same rate while rate of demethylation reduces due to the decrease in active methylated receptor 𝑋 K F caused by change in attractant concentration. This cause increasing the number of methylated receptors including the growth in the number of active receptors. The rate of change of methylated receptors 𝑋 K C +𝑋 K F after the change remain the same. Consequently, the new steady state activity level stays the same. Therefore, the activity after change in attractant concentration is equal to the activity before changing attractant concentration. i(j J k l j J m ) in = 𝑉 Z R− o p ] j J m ql j J m → 𝐴 e = 𝑋 K F = qo r Z o p ] Ho r Z → 𝐴 e = 𝐴 ` (2.17) The bacterium adapts itself to the attractant concentration changes. The reason is that the total number of active methylated receptors adjusts itself due to the balance in demethylation rate. This means, independent of the concentration of attractant, the total number of active methylated receptor return to its fixed point 13 [98], [99]. In this model, exact adaptation happens for a wide range of parameters while the value of steady state activity, which the bacterium adapts itself, depends on the parameters. This mean the exact adaptation is a robust feature of the model, while the steady-state activity level is a fine tuned feature [67], [84], [100], [101]. Figure 2.4 Two simplified models for methylation and demethylation process in E. coli bacteria, which capture exact adaptation property. a) Fine-Tuned model. In this model, the methylated receptors are active; while the unmethylated receptors are inactive. Therefore, exact adaptation depends on proportion of different biochemical parameters. b) Barkai-Leibler model. In this model, unmethylated receptors are in inactive state; while the methylated receptors can transit between active and inactive states. Therefore, exact adaptation is robust and it is valid for a wide range of biochemical parameters 2.5 Combination of The Mentioned Models to Simulate Standard Model of E. coli. In this section, we provide a more detailed description of the full chemical pathway model of E. coli used in our BNSim simulation environment. In general, the chemotaxis pathway of E. coli has three components: the cooperative chemoreceptor, the phosphorylation pathway, and the flagellar motor [83], [102]–[105]. In what follows, we describe these three components of the chemotaxis pathway. 2.5.1 Cooperative Chemoreceptor Cluster We use a combination of introduced models for cooperativity and methylation kinetics to simulate the chemoreceptors of E. coli. To describe the cooperativity of chemoreceptors, we use the latest MWC model that captures the essential features of bacteria chemotactic activity and adaptation [104]. More specifically, each functional methyl-accepting chemotaxis protein (MCP) receptor complex can be either in the active or the inactive states, which are determined by a free energy difference F(m, [L]), where m is the methylation level of the receptors, [L] is the concentration of the ligand, and N is the number of receptor dimers in a MCP complex. We use N = 6 for Tar receptors (for aspartate) in a complex. To capture the 14 essential methylation and demethylation noise, we model 370 receptor clusters individually [106]. The average activity of the receptor can be expressed as [66], [107]: 𝑎 = (1+exp 𝑁×𝐹 𝑚, 𝐿 ) Ha (2.18) According to the MWC model, the free energy difference can be written as: 𝐹 𝑚, 𝐿 = 𝑓 K 𝑚 −𝑙𝑛 1+ 𝐿 𝐾 F +𝑙𝑛 1+ 𝐿 𝐾 C (2.19) In equation (2.17), 𝑓 K (𝑚) is the methylation level dependent free energy difference, and 𝐾 F and 𝐾 C are the dissociation constants of the ligand to the active and the inactive receptors, respectively. For the binding of L-aspartate onto the Tar receptors, we use the values fitted to in vivo FRET (fluorescence resonance energy transfer) data [66], namely, 𝐾 F = 3 mM and 𝐾 C = 18.2 𝜇𝑀. The impact of receptor methylation on its free energy is considered to be a linear function of m as suggested by recent experimental work [95], [108]: 𝑓 K 𝑚 = 𝛼 𝑚 ` −𝑚 (2.20) In equation (2.18), 𝛼 = 1.7 and 𝑚 ` = 2. The methylation kinetics is based on the Barkai and Leibler model for a near perfect adaptation system [69]. More precisely, the methylation kinetics is assumed to have a linear form: 𝑑𝑚 𝑑𝑡 = 𝑘 Z 1−𝑎 −𝑘 ] 𝑎 (2.21) In equation (2.21), 𝑘 Z = 𝑘 ] are the methylation and demethylation rates, respectively. Figure 2.5 shows the simulation results of bacterium chemosensitivity and its precise adaptation. Bacteria adaptation is a key feature of the chemotaxis system, which shows its ability to return to the same level of chemoreceptor activity over a wide range of chemical gradients. Figure 2.5 Simulation results for chemoreceptor activity in E. coli showing precise adaptation. The figure is obtained by averaging simulation results over 370 clusters of receptors [106]. 15 2.5.2 Phosphorylation Relay and Flagellar Motor An active receptor enhances the autophosphorylation of the receptor-associated kinase CheA, which transmits the signal to the flagellar motors by the phosphorylation of a diffusive response regulator CheY. In most chemotaxis full pathway models, the concentration of phosphorylated CheY-P is assumed to be proportional to the kinase activity, 𝑌 = 𝛽𝑎 (𝑡), without considering the nonlinear dependency [107]. Thermal fluctuations and upstream signaling cause the flagellar motor to spontaneously changes between CCW and CW states. To model the state transitions, the SPEC model uses a Hill function to calculate the probability of tumble with 𝑌 , and assumes an average fixed tumble time 𝜏 = 0.2 s [66]. Specifically, when a bacterium is running, the probability of the cell going into a tumble state is 𝑝 𝑎 = 𝜏 a Ha ( F F ) , where 𝑎 a/e is a fitted constant, and H is the Hill coefficient of the motor response function. However, under some gradient conditions, the bacteria flagellar motor may stay in the tumble state for a time significantly shorter or longer than the average value, which biases the swimming dynamics considerably [109]. Therefore, we adopt a two state potential well model to describe the motor behavior of E. coli, which sets their two states in two potential wells. The energy barriers of CCW to CW and CW to CCW transitions are 𝐺 ` ( 𝑌 ) and -𝐺 ` ( 𝑌 ), with transition rates 𝑘 H and 𝑘 l , respectively [102]: 𝐺 ` 𝑌 = 𝑔 ` 4 − 𝑔 a 2 ( 𝑌 𝐾 + 𝑌 ) (2.22) 𝑘 l = 𝜔 ` exp (𝐺 ` 𝑌 ) (2.23) 𝑘 H = 𝜔 ` exp (−𝐺 ` 𝑌 ) (2.24) Where parameters 𝜔 ` = 1.3 𝑠 Ha , 𝑔 ` = 𝑔 a = 40 𝑘 ] 𝑇, and 𝐾 = 3.06 𝜇𝑀 are chosen to fit the experimental data. 2.5.3 Bacteria Transport Model: In a bounded channel with chemical gradient, change in the motion of single bacteria leads to a biased distribution of bacterial density. Based on probabilistic modeling of single bacteria, the model described in this section relates the population kinetics to the motility parameters of individual cells. The chemotactic velocity (𝑉 ) introduced in the model quantifies the effect of chemotaxis on the bacterial population transport; by construction, 𝑉 is directly related to the tumble rate bias of individual bacteria. Therefore, by extracting 𝑉 from this model and determining it through tracking of single bacteria, we quantify chemotaxis from the measurement of individual cells, as opposed to most recent studies, which have made measurements at the population level [66], [110]–[112]. 16 In an environment with a one-dimensional (along X-direction) chemoattractant gradient, the transport kinetics of bacteria density can be described with the following equation [113]: 𝐽 = −𝜇 𝜕𝐵 𝑥,𝑡 𝜕𝑥 +𝑉 𝐵(𝑥,𝑡) (2.25) In equation (2.25), 𝐽 is the density flux of bacteria through a slice that is perpendicular to the X axis, 𝜇 is the motility coefficient, B(x,t) is the bacteria density, and 𝑉 is the chemotactic velocity as previously discussed. The motility coefficient 𝜇, similar to the diffusion coefficient, is a proportional measurement between the bacterial density flux and the gradient of the bacterial density. Bacteria transport in the channel is subject to two effects: the diffusivity motility and the chemotactic drift, which are described by the first and second terms on the right hand side of equation (2.25), respectively. The diffusivity motility keeps the transport bacteria down the gradient of bacterial population density, while the chemotactic velocity pumps bacteria up the chemoattractant gradient. From the probabilistic modeling of individual bacteria, the motility coefficient and chemotactic velocity can be expressed as follows [104], [114]: 𝑉 = e( H ) ( l ) 𝜐 (2.26) 𝜇 = 2𝜐 e 3(𝑝 H +𝑝 l )(1−𝑐𝑜𝑠 < 𝜙 >) (2.27) In equation (2.26) and (2.27), 𝑝 l and 𝑝 H are the average tumble rates when bacteria travel up and down the chemical gradient, respectively. 𝜐 is the average 3D swimming speed of the bacteria. < 𝜙 > is the average tumble angle, which is defined as the swimming direction change during a tumble. Instead of a simple model based on average rate processes, models based on random walk theories that consider distributions of running and tumbling durations are also available [115], [116]. These advanced models could provide a more accurate characterization of single bacteria motility and population transport if both the run and tumble durations of E. coli can be measured experimentally in the future. 17 3 Mathematical Characterization of Single Bacterium Motion 3.1 Introduction In this chapter, we investigate single bacterium motion characteristics while swimming in a three- dimensional environment (the content in this chapter have been published in [117]). We take into account realistic conditions like volume exclusion (i.e., no two bacteria can occupy the same space at the same time), chemical interactions among bacteria 3 (e.g. autochemotaxis), obstacles [118], [119], and heterogeneous distribution of chemoattractant in the environment and investigate non-ergodicity, nonlinearity and multi-fractality of single bacterium motion. We identify fractal characteristic of a single bacterium motion, which could have a fundamental impact on mathematical modeling of both single bacterium and population swimming dynamics. We analyze the motion trajectories of E. coli and S. marcescens; these are two types of bacteria which have very similar swimming dynamics, behavior, and chemotaxis [20], [21]. Extracting single bacterium trajectories of motion from computer simulations for E. coli and S. marcescens (Figure 3.1a-b show the different simulation setups that we used for the chemoattractant gradient in the environment) and applying single-cell visual tracking methods to real experiments on S. marcescens (Figure 3.1c) reveal the fractal characteristic of bacteria motion. To elucidate the complexity of the bacterial dynamics, we investigate the individual trajectories of simulated E. coli and S. marcescens beside trajectories of real S. marcescens at multiple scales in space and time; this way, we show single bacterium motion has a phase transition from a superdiffusive behavior to a normal diffusion and lastly to a subdiffusive pattern based on bacteria density and chemoattractant distribution. We also show that bacterium motion exhibits a multi-fractal behavior. The multi-fractal feature means that bacterial motion is self-similar, which is a concept dating back to the pioneering work of Kolmogorov to explain the chaotic behavior [22], [23]. Consequently, the multi-fractal formalism offers mathematical techniques to capture this nonlinear self-similarity in their motion. The remainder of this chapter is organized as follows: In sections 3.2 and 3.3, we briefly overview the simulation setup and experimental setup to generate the bacterium motion trajectories. Next, we investigate 3 Autochemotaxis happens when bacteria excrete the converted substrate succinate molecules in the environment into chemoattractant aspartate molecules [109]. 18 the swimming dynamics of simulated E. coli. We explore the second order moment of motion using an ergodicity test under different conditions including bacteria density and chemoattractant distribution. Then we explore the third order moment of motion with a nonlinearity test. Next, we discuss that focusing on second and third order moment of a single bacterium’s motion cannot capture the entire spatial and temporal complexity of the motion; therefore, we perform multi-fractal analysis which considers higher order moments of their motion. We show that the simulated E. coli bacterial swimming displays multi-fractal characteristics. Then, in section 3.5, we perform the same analysis for swimming dynamics of simulated S. marcescens. In section 3.6, using real experiments on S. marcescens in vitro, we show that their motion is multi-fractal, validating our conclusions for simulated bacteria. 3.2 Simulations For completeness of our analysis, we used computer simulations of known mathematical models for the chemical pathways of bacteria E. coli and S. marcescens to capture their dynamics through simulations in addition to experimental data. The main reasons for using simulations beside real experiments is the ergodicity test. For performing ergodicity test and measuring the time averaged mean square displacement and ensemble mean square displacement, we need at least 1000 trajectories of exactly the same bacteria starting with exactly the same initial condition. Due to constraints on experimental conditions, we are not able to perform the same experiment 1000 times exactly with the same initial conditions to track exactly the same bacteria in all of them. To tackle this issue, we used BNSim [106], an open-source, parallel multi- scale stochastic modeling platform integrating various simulation algorithms with genetic circuits and chemotaxis pathway modeled in a complex three-dimensional environment (More details have been explained in chapter 0 (section 2.5), the source code can be downloaded from http://www.ece.cmu.edu/~sld/bnsim/index.html). Specifically, to simulate the chemoreceptors, we used the Monod-Wyman-Changeux model in which the receptor homo-dimers assemble into fully cooperative signaling teams that switch rapidly between active and inactive states [66]. The methylation kinetics is based on the celebrated Barkai and Leibler model for a near perfect adaptation system [95], [120]. For the signal transduction from chemoreceptor to the flagella motor regulator Yp, the concentration of phosphorylated CheYp is assumed to be proportional to the kinase activity without considering the nonlinear dependence [121]. Finally, we used a two-state model to describe the motor behavior of bacteria, which sets the clockwise (CW) and counterclockwise (CCW) states in two potential wells [122], [123]. The transition rates have been fitted to experimental data of E. coli [103] and S. marcescens [20], [21]. We considered autochemotaxis in all of our simulations. Autochemotaxis happens when bacteria excrete the converted substrate succinate molecules in the environment into chemoattractant aspartate molecules [109], [124]–[126]. Once the succinate gets depleted, bacteria start consuming the chemoattractant aspartate 19 excreted by themselves and they may detect each other’s presence through this process [15]. Of note, this form of chemical interaction has been defined as cue in the chapter entitled “communication in bacteria’ by Diggle and his coworkers in [127]. The simulations have been done for three different cases: Case 1. To simulate E. coli, we used a cubic environment with dimension 5000x5000x5000 µm 3 (Figure 3.1a). The whole environment was subdivided into 10 6 smaller cubes each with size of 50x50x50 µm 3 . Simulated bacteria had a constant run speed of 14 µm/s and a rotational diffusion constant of 0.062 rad e .s [93], [128]. After each tumble, cells were reoriented in a new direction that is randomly sampled from a gamma distribution with scale parameter 18.32 and shape parameter 4 matching the experimentally observed distribution of new run angle. In this set of simulations, different population densities of 8x10 2 bacteria/cm 3 , 8x10 3 bacteria/cm 3 , and 8x10 4 bacteria/cm 3 in the environment were considered. These population densities correspond to total number of 10 2 bacteria, 10 3 bacteria, and 10 4 bacteria respectively. To study the effect of chemotaxis in the environment, we considered a steady state linear gradient (10 -4 mM/µm) of chemoattractant (aspartate) in y-direction. We simulated each bacteria density considering two different condition: with the gradient of chemoattractant in the environment and without the gradient of chemoattractant in the environment. To analyze the effect of obstacles on bacteria motion, we considered 0.01 percent of the small cubes uniformly distributed in the environment to be impenetrable (i.e. totally 10 4 number of cubes are impenetrable), thus when the bacterium hits one of these cubes it cannot continue to enter the cube and it will continue tumbling to choose a new direction for run. Case 2. In this case we had the same set up for environment and simulated bacteria as case 1 to simulated E. coli. The only difference was in the population densities we simulated and chemoattractant gradient in the environment. In this case, we considered a wider range of population densities than case 1, namely, 8x10 2 bacteria/cm 3 , 8x10 3 bacteria/cm 3 , 8x10 4 bacteria/cm 3 , 8x10 5 bacteria/cm 3 , and 8x10 6 bacteria/cm 3 in the environment. These population densities correspond to total number of bacteria from 10 2 bacteria up to 10 6 bacteria in the simulations. The other significant difference between this case and case 1 was that we considered a steady state radial gradient (10 -4 mM/µm) of chemoattractant (aspartate) from the injection location towards the targeted region, namely, ( 𝑥 e +𝑦 e +𝑧 e `.¡ ×10 H¢ mM/µm). We simulated each bacteria density considering two different condition: with the gradient of chemoattractant in the environment and without the gradient of chemoattractant in the environment. Case 3. To simulate S. marcescens, we used a three-dimensional environment with dimension 10000x500x150 µm 3 . The whole environment was subdivided into smaller cubes each with size of 50x50x50 µm 3 . For the simulation, densities of 10 4 bacteria/cm 3 , 10 5 bacteria/cm 3 , 10 6 bacteria/cm 3 , and 10 7 20 bacteria/cm3 in the environment were considered. These population densities correspond to total number of 7 bacteria, 75 bacteria, 750 bacteria, and 7500 bacteria respectively. To study the effect of chemotaxis in the environment, we considered a steady state linear gradient (10 -4 mM/µm) of chemoattractant (aspartate) in the y-direction. We simulated each bacteria density considering two different condition: with the gradient of chemoattractant in the environment and without the gradient of chemoattractant in the environment. Figure 3.1 Simulation configuration of E. coli in BNSim (see Methods) and single particle tracking of S. marcescens in vitro experiments. (a) In case 1 for simulation of E. coli in BNSim, a cubic environment with dimensions 5000x5000x5000 µm 3 has been considered. The whole environment subdivided into 10 6 smaller cubes each with size of 50x50x50 µm 3 . We consider a linear gradient (10 -4 mM/µm) of chemoattractant (L-aspartate) only in the y-direction and no gradient in other directions. (b) In case 2 for simulation of E. coli in BNSim, we consider the same configuration as case 1 for bacteria and environment. We only changed the gradient of chemoattractant in the environment. In this case we consider a linear gradient of chemoattractant (L-aspartate) from the injection location to the targets, meaning there is linear gradient in 3 different dimensions. (c) Single particle tracking of S. marcescens in an in vitro environment with dimensions of 10000x500x150 µm 3 with linear gradient (10 -4 mM/µm) of chemoattractant (L-aspartate) in the second direction (y-direction). 21 3.3 Experiments Laboratory experiments were conducted to measure the swimming behavior of live and motile S. marcescens under two conditions: (1) an isotropic environment (i.e. uniform distribution of substances) with dimension equal to 10000×500×150 µm and (2) a stable-linear gradient of L-aspartate, which is a canonical chemoattractant for bacterial species such as S. marcescens and E. coli. 3.3.1 Bacteria Culture To conduct the experiments, S. marcescens (ATCC 274, American Type Culture Collection, Manassas, VA) was cultured using established protocols [21], [129]–[131]. The bacteria were first grown to the exponential state in a liquid culture (25 g Difco LB Miller Broth and 1 L deionized (DI) water, pH 7.0) on a shaker at 37°C for 3.5 - 4 hours. Then an aliquot of 2.0 µL of the liquid culture was placed on an agar plate (25 g Difco LB Miller Broth, 6 g Bacto Agar, 5 g glucose, 1 L de-ionized water), and the agar plate was incubated at 30°C for 16 - 20 hours. For the experiments, bacteria from the leading edge of the colony were extracted and then diluted in motility buffer (0.01 M KH 2 PO 4 , 0.067 M NaCl, 0.1 mM EDTA, pH = 7.0) to appropriate densities for the measurements. Subculturing the bacterial cells on an agar plate allows the most motile bacterial cells, which have been shown to be located along the leading edge of a spreading colony [21], to be selectively chosen for the study. 3.3.2 Experimental Setup A three-channel microfluidic concentration gradient generator was fabricated to measure bacterial swimming parameters under chemotaxis and without chemotaxis; detailed descriptions of the device can be found in [20]. The gradient generator consisted of three parallel microfluidic channels within an agarose gel. By controlling the concentration of the chemoattractant in the outer two channels of the device, a linear chemoattractant concentration profile could be generated by molecular diffusion in the middle channel in which the bacterial suspension was placed. The gradient generator was calibrated by observing the diffusion of 10 -4 M fluorescein (Sigma-Aldrich Co.) in the device as described in [20] and [21]. Fluorescent intensity profiles from these studies revealed that a stable, linear gradient developed within 20 minutes. Given that the diffusion coefficient of L-aspartate is 9.0 × 10 -6 cm 2 /s and the diffusion coefficient of fluorescein is about half that of L-aspartate [20], it can be concluded that about 10 minutes is required to develop a stable chemoattractant concentration gradient of L-aspartate and the gradient remains stable beyond 10 minutes. In our study, the bacterial swimming motion was measured after 15 minutes, ensuring a stable linear concentration gradient within a quiescent fluid environment. 22 3.3.3 Experimental Conditions Chemotaxis of S. marcescens was measured under an L-aspartate gradient of 0.2 mM/mm, with an average concentration of 0.1 mM. To measure the bacterial swimming behavior in an isotropic environment, we used the same device but with a uniform distribution of motility buffer (without an L-aspartate gradient). For both cases, with and without the chemoattractant gradient, the bacterial swimming parameters were studied for three different cell densities, which were 10 6 bacteria/cm 3 , 10 7 bacteria/cm 3 and 10 8 bacteria/cm 3 . All the experiments were conducted at room temperature (19 – 22°C). 3.3.4 Imaging and Tracking The bacteria were imaged with an inverted microscope (40x objective, Zeiss AxioObserver.A1, Carl Zeiss, Oberkochen, Germany), and videos were obtained at 88 frames per second with a CCD camera (FOculus FO134SB). All video data was acquired at least 20 µm (~ 20 body lengths) from any surface to minimize the wall effects on the bacterial swimming behavior, and no obstacles hindered the movement of the bacteria (except for other bacterial cells). The bacteria were tracked in three-dimensions by an in-house visual tracking program developed in Matlab (R2012a, The MathWorks, Inc, Natick, MA). The tracking program extracts the x-y positions of bacteria in the video frames using image thresholding techniques, while the z- position was obtained by making use of the lens aberrations created by out-of-focus bacterial cells. This defocused optical tracking method has been applied by several groups to track the three-dimensional motion of bacteria and small particles [21], [132]–[135]. A linear relationship exists between the radius of the aberration ring appearing around an out-of-focus cell and the cell’s vertical distance (along z-axis) to the focal plane. Equation (3.1) exhibits this relation. 𝑟 = 𝑐 a 𝑧 +𝑐 ` (3.1) In equation (3.1), c 1 and c 2 were calibrated by performing an experiment in which the bacteria were fixed a known distance away from the focal plane. The variables were determined from [21]. Therefore, the z- position of a bacterial cell could be determined by measuring the size of the aberration ring around the cell. To connect bacteria positions in individual frames and form the time trajectories of bacteria up to tens of seconds in length, we used a nearest neighbor algorithm. In other words, we connected the disjoint path segments by using a weighting function, which determined the likelihood that two path segments were produced serially by the same bacteria. Equation (3.2) explains this weighting function in more detail. This algorithm helps to connect positions of bacteria in different frames. 𝑤 = 𝑐 e 𝑑 e +𝑐 ∆𝜈 e +𝑐 ¢ ∆𝜃 +𝑐 ¡ ∆𝑡 (3.2) In equation (3.2), ∆𝜈 is the change in speed between the end of the first trajectory and the start of the second trajectory, ∆𝜃 is the change in heading during a tumble, and ∆𝑡 is the change in time. d is the difference in 23 position between the start of the second trajectory and where as object moving at the final speed and heading of the first trajectory would be at a time ∆𝑡 after the end of that trajectory. The coefficients c i are empirically fitted constants. Further details about the tracking algorithm are described in [21]. 3.4 Investigation of Simulated E. coli Swimming Dynamics 3.4.1 Introduction on Ergodicity Test Simply speaking, if the dynamics of a bacterium’s motion is nonergodic, the time averaged mean square displacement (TAMSD) (equation (3.3)) is not equal to the ensemble mean square displacement (EMSD) (equation(3.4)). In mathematical terms, the system is nonergodic if TAMSD ≠ EMSD [136]–[138]. In the TAMSD formula, 𝛥 and 𝑇 represent the lag time and the overall measurement time, respectively. In the EMSD formula, 𝑃(𝑟,𝑡) is the probability density function to find bacteria at position r at time t (in TAMSD formula the overline denotes the time average). lım →® δ e (Δ = t,T) = 1 T−Δ r (t+Δ) −r (t) e HS ` dt (3.3) < r e (t)> = r e P r,t d r= < |r (t) −r (0)| e > (3.4) 3.4.2 Exploring the Second Order Moment of Bacterium Motion with the Ergodicity Test The dynamics of bacterial population is affected by numerous factors such as the spatio-temporal gradient of chemoattractant within the environment [139], [140], bacteria density, volume exclusion [141], as well as the intensity of autochemotaxis [142], and saturation of bacteria receptors. One important consequence of the mentioned dependencies is that the collective motion of bacteria is nonergodic in nature [143]–[145]. To investigate this, we performed an ergodicity-breaking test which checks that long-time averages differ from ensemble averages of bacterium displacement. Figure 3.2a-d summarize the ergodicity investigations on the simulated E. coli swimming dynamics under various bacterial densities, chemoattractant gradient, and autochemotaxis effects. For instance, Figure 3.2a-b show the difference between TAMSD and EMSD plots for a simulated E. coli bacterium in a population of 8x10 4 bacteria/cm 3 in an environment with a linear gradient of chemoattractant and without obstacles except other bacteria. More precisely, the TAMSD for this case when the lag time 𝑡 = 10 (s) is 𝛿 e 𝑡 = 10,𝑇 = 1000 = 15631 (µm e ), while EMSD obtained over 100 trajectories for the same timing conditions is < 𝑟 e (𝑡 = 10 (s),𝑇 = 1000(s))> = 12961 (µm e ). The difference between TAMSD and EMSD in this case is 20 percent. The TAMSD for a lag time 𝑡 = 10 (s) is independent of initial time reference, which means that bacteria will move a distance equal to 125.02 (µm) on average. In contrast, the 24 EMSD for the same time t = 10 (s) implies that bacteria will travel (on average) a distance of 113.84 (𝜇𝑚) which has a 10 percent difference. The significant difference between the TAMSD and EMSD shows that for simulated E. coli, the motion is nonergodic and time dependent [145]. Comparing the slopes of the lines fitted to both TAMSD and EMSD plots confirms that simulated E. coli motion is nonergodic (the slope of the line fitted to the TAMSD plot is 1.198 (+/- 0.007) while that for the EMSD plot is equal to 1.290 (+/- 0.009)). For completeness, we test all cases for different points in time (Table 3.1). The cases we analyzed cover different bacteria density and chemoattractant distribution in environment. In all of these cases the equality between TAMSD and EMSD does not hold; this confirms the nonergodicity of the bacteria motion. In other words, this nonergodicity holds independent of bacteria population dynamics and chemoattractant distribution in the environment. Table 3.1 Comparing TAMSD and EMSD for different densities of cases 1 for simulated E. coli bacterium trajectory Case 1 Time lag TAMSD EMSD Error 8x10 2 bacteria/cm 3 , with chemoattractant 10 (s) 16529 13261 % 24 8x10 2 bacteria/cm 3 , with chemoattractant 100 (s) 736230 719880 % 2 8x10 2 bacteria/cm 3 , with chemoattractant 400 (s) 3779700 3499200 % 8 8x10 2 bacteria/cm 3 , with chemoattractant 800 (s) 8796900 7886600 % 11 8x10 3 bacteria/cm 3 , with chemoattractant 10 (s) 15787 19082 % 20 8x10 3 bacteria/cm 3 , with chemoattractant 100 (s) 741230 664580 % 11 8x10 3 bacteria/cm 3 , with chemoattractant 400 (s) 2420200 3220700 % 33 8x10 3 bacteria/cm 3 , with chemoattractant 800 (s) 11180000 6648500 % 40 8x10 4 bacteria/cm 3 , with chemoattractant 10 (s) 15631 12961 % 20 8x10 4 bacteria/cm 3 , with chemoattractant 100 (s) 597950 687400 % 14 8x10 4 bacteria/cm 3 , with chemoattractant 400 (s) 3706100 3395200 % 9 8x10 4 bacteria/cm 3 , with chemoattractant 800 (s) 3979800 5281200 % 32 25 Figure 3.2 Mean square displacement (MSD) plots for simulated E. coli in BNSim. (a) Time averaged MSD of one E. coli bacterium motion simulated in BNSim as a function of lag time Δ (δ e (Δ,T) = a HS r (t+Δ) −r (t) e HS ` dt) for a density of 8x10 4 (bacteria/cm 3 ) without chemoattractant and obstacles. This plot shows superdiffusion behavior (MSD (t) = k ± t ± ,1< α< 2). (b) Ensemble averaged MSD of one E. coli motion simulated in BNSim vs. time t calculated from the following formula (MSD (t) = < |r (t) −r (0)| e >) over 100 trajectories for a density of 8x10 4 (bacteria/cm 3 ) without chemoattractant and obstacles. This plot also shows superdiffusion. The observed differences between these two plots (a and b) show an ergodicity breaking in the system and demonstrate the nonergodicity of bacteria dynamics. (c) The TAMSD of the center of whole swarm for E. coli motion simulated in BNSim shows the effect of population density and chemoattractant on center of swarm motion (case 2). (d) The TAMSD of one E. coli trajectory simulated in BNSim for different cases show the effect of chemoattractant and density on a single bacterium motion (case 2). By increasing the population density, we observe a phase transient in bacteria behavior from superdiffusion to normal diffusion and then subdiffusion (case 2). The nonergodic behavior of the motion can be also recognized from the anomalous diffusivity of EMSD plots. Generally speaking, the anomalous diffusion is a diffusion process characterized by a nonlinear relationship of the EMSD with time [146]–[150]. In contrast, a normal diffusion process has an EMSD that varies linearly with time. Figure 3.2b shows bacteria driven motion displays a superdiffusion type of behavior [151], which can be modeled using a power law form of the second moment with exponent ranging between 1 and 2 (𝑖.𝑒.< 𝛿 e (𝛥,𝑇)>= 𝐾 ³ 𝛥 ³ ,1 ≤ 𝛼 ≤ 2) for different bacteria density cases for the simulated E. coli. 26 Figure 3.2c-d show that single bacterium trajectory and the movement of the center of the group exhibit superdiffusion behavior for some cases in our simulations. The superdiffusion behavior predicts that its corresponding random walk may be characterized by a Levy distribution [152], [153]. To further elucidate this behavior of simulated E. coli, we investigate the histogram of the distance traveled by bacteria 𝑃(𝑟,𝜏) in a specific lag time 𝜏 (Figure 3.3a) and fit two type of Gaussian and Levy type distribution to it. We observe that it exhibits a Levy walk type of distribution and does not obey Gaussian distribution. The origin of this behavior is the lack of nutrient concentration in our simulation model, which leads to a high level of noise in the methylation model. This noisy fluctuation in methylation model is the reason for the power- law run length distribution in simulated bacterium motion, which results in a transition from classical random walk to Levy walks in single bacterium motion [152]. Figure 3.3a shows that the side peak of the bacterium displacement histogram 𝑃(𝑟,𝜏) shifts towards right hand side and broadens with time. Figure 3.2c-d show the effect of population density on TAMSD for different bacteria density cases for simulated E. coli. The TAMSD variation of the center of the entire bacteria population as a function of time lag 𝛥 (Figure 3.2c) shows that increasing the bacterial density makes the population less motile. Volume exclusion effects and chemical interactions caused by autochemotaxis between bacteria are the main reasons of this behavior. Figure 3.2d shows the TAMSD of a single bacterium, which is in agreement with the results obtained from the TAMSD of the entire swimming trajectory. One reason for this behavior is that increasing their population causes a significant increase of the autochemotaxis. In turn, this may cause receptor saturation which can “confuse” bacteria in choosing its preferred direction of motion in the environment; this means that their motion may become restricted due to their high accumulation in the environment. Hence, TAMSD of bacteria motion decreases (for the same time interval) with the increase of bacteria density. Figure 3.2d magnifies the volume exclusion effect for the case of 8×10 µ bacteria/cm 3 (purple and magenta lines). As one can see, in this case, the bacterium is not able to move and therefore it remains in the neighborhood of its original location. By increasing the bacteria population up to 8×10 ¡ bacteria/cm 3 , we observe a phase transition in their motion from a superdiffusive behavior to a normal diffusion and lastly to a subdiffusive pattern (Table 3.2). Table 3.2 Slope of line fitted to TAMSD plot of a single bacterium Bacteria density Chemoattractant condition 𝛼(±𝑒) Diffusion type 8x10 2 (bacteria/cm 3 ) Without chemoattractant 1.402(±0.026) superdiffusion 8x10 2 (bacteria/cm 3 ) With chemoattractant 1.483(±0.022) superdiffusion 8x10 3 (bacteria/cm 3 ) Without chemoattractant 1.229(±0.015) superdiffusion 8x10 3 (bacteria/cm 3 ) With chemoattractant 1.220(±0.022) superdiffusion 8x10 5 (bacteria/cm 3 ) Without chemoattractant 0.79(±0.020) subdiffusion 8x10 5 (bacteria/cm 3 ) With chemoattractant 1.07(±0.010) superdiffusion 8x10 6 (bacteria/cm 3 ) Without chemoattractant 0.040(±0.010) subdiffusion 8x10 6 (bacteria/cm 3 ) With chemoattractant 0.050(±0.010) subdiffusion 27 Besides the TAMSD plots, the displacement histogram 𝑃(𝑟,𝜏) shows a higher probability for a longer displacement in the same time interval related to lower bacterial densities (Figure 3.3b). When constraints like other bacteria in the environment limit a bacterium’s motion, then the bacterium tends to oscillate more and randomly move in the environment rather than moving directionally. Figure 3.3 Effects of time, bacteria population density and chemoattractant on E. coli motion simulated in BNSim (a) Effect of time on bacteria motion. Log-log plot of P(r,τ) as a function of distance r for different lag times τ = 100 s,500 s,1000 s. This result is for 8x10 2 bacteria/cm 3 in an environment without chemoattractant and obstacles. This plot shows that the bacteria move more with time and get further away from their initial locations. The peak of the plot decreases more (case 1). (b) Effect of density on bacteria motion. Log-log plot of P(r,τ) for lag time τ = 1000 s compared for different bacteria (E. coli) density simulated in BNSim in an environment with chemoattractant gradients and without obstacles. By increasing the bacteria population, their motion will be restricted and they are able to move less freely in the environment (case 1). (c) Effect of chemoattractant on the bacteria motion. Log-log plot of P(r,τ) for lag time τ = 1000 s compared for the case with and without chemoattractant in the environment for the population of 8x10 2 bacteria/cm 3 . By adding chemoattractant in the environment the bacteria tend to move more and oscillate less so they get further away from their initial conditions in the direction of increasing gradient of chemoattractant (case 1). 28 We also study the effect of L-aspartate gradient on TAMSD of bacteria motion trajectory. Figure 3.2c-d show that by adding L-aspartate gradient (case 2) to the environment, independent of bacteria density, the TAMSD will increase for both a single bacterium and entire population for simulated E. coli. Moreover, Figure 3.3c presents the histogram of displacement of bacteria after lag time 𝜏 (i.e. 𝑃(𝑟,𝜏)) for the case with/without L-aspartate gradient in the environment. This figure shows that the probability of covering longer displacements (after the same time lag) by bacteria is higher in the case of L-aspartate gradient in the environment, meaning that bacteria tend to move directionally along a targeted direction and oscillate less when they sense a chemoattractant gradient in their environment while their instantaneous speed stays the same (Figure 3.3c). Up to this point, we focus on nonergodicity of single E. coli motion by investigating the second order moment of their motion under different environmental condition. In what follows, we shift our focus to the third order moment of their motion by performing a nonlinearity test. 3.4.3 Introduction on Nonlinearity Tests 3.4.3.1 Mardia’s Test To check the nonlinearity of a bacterium motion we performed Mardia's test on the data set for bacteria trajectory. This method is based on multivariate extension of skewness measure. The bacterium motion is nonlinear if the skewness of trajectory of its motion is non-zero. This also demonstrates bacterium motion does not have a normal distribution [154]–[156]. Based on Mardia’s approach, the skewness of the data for a sample of {r 1 ,r 2 ,... r n } of k dimensional vector relates to the third moment of the data set. The first step to find the third moment of data set is calculating the covariance matrix 𝛴 from equation (3.5). 𝛴 = 1 𝑛 (𝑟 ¹ −𝑟) (𝑟 ¹ −𝑟) º » ¹¼a (3.5) In equation (3.5), 𝑟 is the mean of the data set. In the next step we calculated skewness from equation (3.6) 𝑠𝑘𝑒𝑤𝑛𝑒𝑠𝑠 = 1 6𝑛 [ (𝑟 C −𝑟) º 𝛴 Ha (𝑟 ¹ −𝑟)] » ¹¼a » C¼a (3.6) 3.4.3.2 Henze-Zirkler's Multivariate Normality Test Henze-Zirkler's Multivariate Normality Test is another approach to test the nonlinearity of data set. This test is based on a nonnegative functional distance that measures the distance between two distribution functions. We calculate this distance from equation (3.7) 29 𝐷 ¾ 𝑃,𝑄 = ∫ Z |𝑃 𝑡 −𝑄 𝑡 | e ×𝑓 ¾ 𝑡 𝑑𝑡 (3.7) P(t) is the characteristic function of the multivariable normality distribution and Q(t) is the empirical characteristic function of data which we want to test. F b (t) is the weight or kernel function calculated from equation (3.8). 𝑓 ¾ 𝑡 = (2×𝜋×𝑏 e ) Ã ×𝑒 |Ä| ×Å (3.8) p is the number of variables and 𝑡 = (𝑡 d 𝑡) `.¡ . The smoothing parameter b is dependent of sample size n and can be calculated from equation (3.9). 𝑏 𝑛 = 1 2 ( 2𝑝+1 4 ) a l ¢ ×𝑛 a l ¢ (3.9) The Henze-Zirkler statistic has a lognormal distribution which is used to compute the null hypothesis probability (equation (3.10)). 𝑊 »,¾ = 1 2 𝑒 H ¾ e × ÇÈ −2(1+𝑏 e ) H e » É¼a × 1 𝑛 𝑒 H ¾ e× al ¾ × Ç » ¹¼a +(1+2𝑏 e ) H e » ¹¼a 𝐷 ¹É = (𝑋 ¹ −𝑋 É )×𝑖𝑛𝑣(𝑆)×(𝑋 ¹ −𝑋 É )′ 𝐷 ¹ = (𝑋−𝑀𝑋)×𝑖𝑛𝑣(𝑆)×(𝑋−𝑀𝑋)′ (3.10) W n,b is weighted L 2 distance, X is data matrix, MX is sample mean vector and S is the covariance matrix normalized by n. Henze-Zirkler statistic test know in literature as T n,b is as fallow, where 1{.} stands for the indicator function: HZ = n × (4 1{S is singular} + W n,b 1{S is nonsingular}) For a multivariate normal data, the test statistic HZ is approximately lognormal distributed. Then this test calculates the mean, variance and smoothness parameter. It will lognormalized the mean and variance and estimates the p-value [157]–[159]. 3.4.3.3 Royston's Multivariate Normality Test To investigate the nonlinearity of a bacterium motion we perform Royston’s marginal method which tests all the p variates for univariate normality with a Shapiro-Wilk statistic, then combines the p dependent tests into one omnibus test statistic for multivariate normality. Royston transforms the p-Shapiro-Wilk statistics into an approximate Chi-squared random variable, with e (1<e <p) degree of freedom. In this method, we estimate the degrees of freedom by taking into account possible correlation structures between the original p test statistics. If the data is multivariate normal, H function in equation (3.11) is approximately Chi- squared distributed. 30 𝐻 = 𝑒 𝑅 ¹ ¹¼a 𝑝 (3.11) Where 𝑅 ¹ = {𝜑 Ha 1 2𝜑 − 1−𝑊 ¹ Ï −𝑚 𝑠 } e 𝑒 = 𝑝 1+ 𝑝−1 ×𝑚𝐶 W j is the related value of the Shapiro_Wilk statistic for the jth variable in a p-variate distribution. g, m and s are calculated from polynomial approximation and 𝜑 Ha and 𝜑 (.) are, respectively, the inverse and standard normal cdf. mC is an estimate of the average correlation among the R j ’s. This Chi-Square distribution is used to obtain the critical or p-value for the MNV test [159]–[166]. 3.4.3.4 Doornik-Hansen Omnibus Multivariate Normality Test Doornik-Hansen Omnibus Multivariate Normality Test is the multivariate version of the univariate omnibus test to check normality of bacterium trajectory based on the transformed skewness and kurtosis [167]–[169]. We used the correlation matrix and the diagonal matrix of reciprocals of the p standard deviations and transformed a multivariate normal into independent standard normal by equation (3.12). 𝑠𝑡 = 𝑍×𝑉× 𝐿 H a e ×𝑉′ (3.12) In this equation, Z is the standard normalized data matrix, V is the eigenvector matrix and L is the eigenvalues diagonal matrix. In this method, skewness is transformed using the D’Agostino procedure and kurtosis is transformed from a gamma distribution to a chi-square. We calculated Doornik-Hansen statistic from equation (3.13), which approximates to a chi-square with 2p degrees of freedom. 𝐷𝐻 = 𝑧 a ×𝑧 a d +𝑧 e ×𝑧 e d (3.13) 3.4.4 Exploring the Third Order Moment of Bacterium Motion with Nonlinearity Tests Details about the first and second moments of the bacteria trajectory (corresponding to mean square displacement of bacteria motion) do not provide enough information regarding the complexity of their dynamics. Consequently, we analyze the higher order moments of a bacterium’s trajectory by investigating the skewness metric. Skewness represents the third order moment of bacterium motion and helps to investigate whether or not the motion increments (i.e., changes in the spatial position of the bacterium in three-dimensional space between two different time instances.) obtained from the bacterium trajectory obey a multivariable normal distribution in the three-dimensional environment. Nonzero skewness shows that 31 probability density function of motion increments does not follow a normal distribution and this is a signature of nonlinearity [156]. Using the Mardia method [154]–[156], we measured the multivariate skewness of our data for simulated E. coli. The nonzero highly variable skewness plots (Figure 3.4) demonstrate that bacteria motion cannot be modeled by a multivariable normal distribution for all the cases under consideration (i.e. with and without chemoattractant and obstacles in the environment); this confirms the nonlinear behavior. To provide a more comprehensive analysis on the existence of nonlinear behavior, we employed several statistical methods like Henze-Zirkler's multivariate normality test [157], [158], [170], Royston's multivariate normality test [159]–[161], and Doornik-Hansen Omnibus multivariate normality test [167], [168]. Table 3.3 summarizes the results for these tests and further demonstrates the nonlinearity of simulated E. coli bacterium dynamics. 32 Figure 3.4 Skewness plots for some of the simulation E. coli trajectories. a) 8x10 2 (bacteria/cm 3 ) in the cubic environment, with/out linear gradient of chemoattractant in y direction, with/out obstacles on the bacteria way. b) 8x10 3 (bacteria/cm 3 ) in the cubic environment, with/out linear gradient of chemoattractant in y direction, with/out obstacles on the bacteria way. c) 8x10 4 (bacteria/cm 3 ) in the cubic environment, with/out linear gradient of chemoattractant in y direction, with/out obstacles on the bacteria way. In all the cases skewness is unequal to zero meaning bacterium displacement does not have normal distribution; therefore, bacterium motion is nonlinear. 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000 Size of data [n] Skewness 8*10 3 (bacteria/cm 3 ) Without chemoattractants, without obstacles 8*10 3 (bacteria/cm 3 ) With chemoattractants, without obstacles 8*10 3 (bacteria/cm 3 ) Without chemoattractants, with obstacles 8*10 3 (bacteria/cm 3 ) With chemoattractants, with obstacles 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000 Size of data [n] Skewness 8*10 2 (bacteria/cm 3 ) Without chemoattractants, without obstacles 8*10 2 (bacteria/cm 3 ) With chemoattractants, without obstacles 8*10 2 (bacteria/cm 3 ) Without chemoattractants, with obstacles 8*10 2 (bacteria/cm 3 ) With chemoattractants, with obstacles a b c 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000 Size of data [n] Skewness 8*10 4 (bacteria/cm 3 ) Without chemoattractants, without obstacles 8*10 4 (bacteria/cm 3 ) With chemoattractants, without obstacles 8*10 4 (bacteria/cm 3 ) Without chemoattractants, with obstacles 8*10 4 (bacteria/cm 3 ) With chemoattractants, with obstacles 33 Table 3.3 Summary of linearity tests on a single bacterium trajectory from simulated E. coli bacteria. These result are based on different tests on the bacterium trajectory from simulated E. coli with density of 8×10 2 (Bacteria/cm3) in a cubic environment with size 5mm×5mm×5mm without chemoattractant in the environment. Number of variables in all the tests are 3 as the bacterium trajectory is in 3D environment, and the sample size is 1001. All the tests results are in agreement that the bacteria motion does not have normal distribution; hence, bacterium has nonlinear motion. Method Variables Result Henze-Zirkler's Multivariate Normality Test Henze-Zirkler lognormal mean: -0.1520616 Henze-Zirkler lognormal variance: 0.1428457 Henze-Zirkler statistic: 45.5201260 P-value associated to the Henze-Zirkler statistic: 0.0000000 With a given significance = 0.050 Data analyzed do not have a normal distribution (P<alpha) Royston's Multivariate Normality Test Royston's statistic: 241.735218 Equivalent degrees of freedom: 1.903974 P-value associated to the Royston's statistic: 0.000000 With a given significance = 0.050 Data analyzed do not have a normal distribution (P<alpha). Doornik-Hansen Omnibus Multivariate (Univariate) Normality Test Asymptotic statistic: 284.3903 P-value associated to the asymptotic statistic: 0.0000 With a given significance = 0.050 Omnibus Doornik-Hansen statistic: 981.4541 P-value associated to the Omnibus Doornik-Hansen statistic: 0.0000 With a given significance = 0.050 Data analyzed do not have a normal distribution (P<alpha). Skewness plot If the data has a normal distribution, the skewness should be equal to zero. Hence, the data is nonlinear. 3.4.5 Introduction on Multi-fractal Detrended Fluctuation Analysis (MFDFA) We performed MFDFA to investigate the higher order moments of a bacterium trajectory of motion to shed light on the underlying structure in bacteria motion [171]. This analysis helps to verify whether the 34 bacterium motion is multi-fractal or mono-fractal. We explain the MFDFA analysis of the time series in the following steps. The first step in MFDFA is checking that the time series has a random walk like structure. If the time series (x) has a noise like structure, it should be converted to a random walk type time series using equation (3.14). 𝑌(𝑖) = [𝑥 É −< 𝑥 >] C É ¼ a , 𝑖 = 1 ,...,𝑁 (3.14) In equation (3.14), <x> is the mean of time series. Next, based on scaling size s, we divided the time series Y to N s non-overlapping subintervals using equation (3.15). This equation returns the integer part of the division. 𝑁 Ó = 𝑖𝑛𝑡 ( 𝑁 𝑠 ) (3.15) Since the length of the time series 𝑁 is not exactly a multiple of the scaling size s, a short part of the time series in the end will remain after dividing the time series Y to N s non-overlapping subintervals. To consider the remaining part of time series in the end and to not disregard this part, we repeated the same division procedure starting from the opposite end of the time series. In other words, we started dividing the time series into N s non-overlapping subintervals form the end point of the time series, and this time a short part of the time series in the beginning of the time series will remain. Therefore, we had 2𝑁 Ó segments overall [172]. The next step is computing the local Root-Mean-Square variation of the time series in each segment using equation (3.16) and (3.17). This quantifies the local fluctuations in the time series. 𝐹 e (𝑠 ,𝜈) = 1 𝑠 {𝑌[(𝜈−1)𝑠 + 𝑖] − 𝑦 Ô (𝑖)} e Ó C ¼ a 𝜈 = 1 ,...,𝑁 Ó (3.16) 𝐹 e 𝑠 ,𝜈 = 1 𝑠 𝑌 𝑁− 𝜈−𝑁 Ó 𝑠 + 𝑖 – 𝑦 Ô 𝑖 e Ó C ¼ a 𝜈 = 𝑁 Ó +1 ,...,2𝑁 Ó (3.17) In equation (3.16) and (3.17), 𝑦 Ô (𝑖) is the polynomial (e.g. linear, quadratic, cubic or higher order polynomial) fitted to each segment 𝜈. In our analysis we considered a linear fit. In the next step, we averaged over all segments to quantify the qth order fluctuation function. In other words, we calculated the mean of q-order RMS for corresponding scaling size s and order q from equation (3.18). 𝐹 Ö (𝑠) = { 1 2𝑁 Ó [𝐹 e (𝑠 ,𝜈)] Ö e eU × Ô ¼ a } a Ö (3.18) 35 We repeated all the steps for several values of the time scales s. Then we determine the scaling behavior of the fluctuation functions by analyzing log-log plots 𝐹 Ö (𝑠) versus 𝑠 for each value of 𝑞. If the time series has a long-range power-law correlation, equation (3.19) will be valid. In other words, mean of q-order RMS (𝐹 Ö (𝑠)) increases as a power-law for large values of 𝑠. 𝐹 Ö (𝑠)≈ 𝑠 (Ö) (3.19) For simplicity we considered that the length N of the time series is an integer multiple of the scale 𝑠, meaning 𝑁 Ó = 𝑁/𝑠: [𝐹 e (𝑠 ,𝜈)] Ö e eU × Ô ¼ a ~ 𝑠 Ö× Ö Ha (3.20) We call the function 𝐻(𝑞) generalized Hurst exponent which captures the fluctuations in the data set. Generally, exponent 𝐻(𝑞) may depend on 𝑞. For Mono-fractal time series, 𝐻(𝑞) is independent of 𝑞, since the scaling behavior of the variance 𝐹 e (𝑠 ,𝜈) is identical for all the segments 𝜈. If small or large fluctuations scale differently, there will be significant dependence of 𝐻(𝑞) on 𝑞, which means the time series is multi- fractal. Another way to distinguish between mono-fractal vs. multi-fractal time series is to first convert 𝐻(𝑞) to the scaling exponent 𝜏 Ö (defined by equation (10)). And then convert 𝜏 Ö to Holder exponent (i.e. singularity strength) h(q) and singularity dimension 𝐷 Ö (defined by equation (3.22) and (3.23)) via a Legendre transform [172]. 𝜏 Ö = 𝑞×𝐻 𝑞 −1 (3.21) ℎ 𝑞 = 𝜏′(𝑞) (3.22) 𝐷 Ö = 𝑞×ℎ 𝑞 −𝜏(𝑞) (3.23) In the case of mono-fractal time series, the scaling exponent 𝜏(𝑞) has a linear dependency with 𝑞. The linear q-dependency of 𝜏(𝑞) leads to a constant ℎ 𝑞 of these time series because ℎ 𝑞 is the tangent slope of 𝜏(𝑞). The constant ℎ 𝑞 reduces the multi-fractal spectrum to a small arc for the mono-fractal time series. In contrast, the multifractal time series has scaling exponents 𝜏(𝑞) with a curved q-dependency and, consequently, a decreasing singularity exponent ℎ 𝑞 . The resulting multi-fractal spectrum is a large arc where the difference between the maximum and minimum ℎ 𝑞 are called the multi-fractal spectrum width. Thus, the width and shape of the multi-fractal spectrum is able to classify a wide range of different scale invariant structures of time series [171]. 36 3.4.6 Exploring Higher Order Moment of Motion Using MFDFA Analysis Traditional analysis of bacteria motion has focused on measuring the root-mean-square fluctuations of the displacement and thus is limited to first and second order moment analysis. However, relying on the second order statistics is not only inconsistent with the observed non-Gaussian behavior, but also it does not capture the entire spatial and temporal complexity structure of the interactions among bacteria and the environment. Consequently, we employ concepts from statistics (e.g., higher order moments) and fractal theory (e.g., multi-fractal spectrum, generalized Hurst exponent) to perform a multi-fractal [173] analysis and quantify the correlation structure and complexity of bacterium trajectories. More precisely, we compute the generalized Hurst exponent as a function of the higher order moment q and the multi-fractal spectrum. The generalized Hurst exponent of a bacterium trajectory measures the scaling and memory properties of the q th order moment of the distribution of fluctuations. For example, the generalized Hurst exponent for q=2 coincides with the autocorrelation function. In addition, the generalized Hurst exponent helps us distinguish between mono-fractal and multi-fractal behavior. If the generalized Hurst exponent proves to be independent of the order q of the moments, then the bacterium motion can be classified as mono-fractal. In contrast, if the generalized Hurst exponent exhibits any dependence on order q, then the bacterium motion is regarded as multi-fractal. Alternatively, we can investigate the existence of a multi-fractal behavior by analyzing the wideness of the multi-fractal spectrum (i.e., W = ℎ (𝑞) ÝÞß −h (q) Ýâã estimating the range of present Holder exponents in the data). A wider multi-fractal spectrum shows that the multi-fractality degree is stronger meaning the structure of the motion is more complex and inhomogeneous. Therefore, a wider range of fractal exponents are needed to cover the complex structure of the bacterium motion. Another measure for multi-fractality degree is the most probable Holder exponent h(q) corresponding to the peak of the multi-fractal spectrum (i.e. h 0 (q)). Low (ℎ ` (𝑞)) shows the motion of bacteria is more correlated and regular in appearance meaning the motion is mono-fractal. Overall, multi-fractal analysis is an applicable tool to characterize the variability and heterogeneity in bacterium motion structure [174], [175]. Figure 3.5a-d present our results from the multi-fractal detrended fluctuation analysis [172], [176] on motion trajectories of simulated E. coli in BNSim for different bacteria density with and without the L- aspartate gradient in the environment. Based on the simulation results (Figure 3.5a), we observe that increasing the density of bacteria decreases the width of the multi-fractal spectrum (𝑊) and decreases (ℎ ` (𝑞)) as two different measures of the degree of multi-fractality. For instance, in the case of chemoattractant in the environment, the wideness (W) of the multi-fractal spectrum is shrinking from 𝑊 = 0.851 to 𝑊 = 0.182 as the population of bacteria increases from 8x10 2 (bacteria/cm 3 ) to 8x10 6 (bacteria/cm 3 ). Similarly, we observe that the most probable Holder exponent h 0 (q) is shifting from 1.06 to 37 0.02 by increasing bacteria density from 8x10 2 (bacteria/cm 3 ) to 8x10 6 (bacteria/cm 3 ). The comparison shows that by increasing bacteria density, the width of the multi-fractal spectrum decreases; this implies that the bacterium dynamics motion transients from a multi-fractal behavior to a mono-fractal behavior. Simply speaking, the bacterium oscillates more than exhibiting a directional motion. This transition translates into a reduction in h 0 (q). The lower h 0 (q) shows that the process is correlated meaning that the underlying structure in bacterium motion becomes more regular by increasing the bacteria population. The increase in bacteria population makes their multi-fractal spectrum asymmetric, which captures the dominance of low- or high-fractal exponents with respect to medium fractal exponents. In the case of 8x10 6 (bacteria/cm 3 ) with chemoattractant in the environment, the right-skewed spectrum shows relatively strong weight for low-fractal exponents. From Figure 3.5a-b, we can conclude that adding chemoattractant to the environment has exactly the opposite effect on a single bacterium’s motion compared to the effect of increasing density. For example, we considered the two cases of bacteria population equal to 8x10 2 (bacteria/cm 3 ) without and with chemoattractant in the environment. The wideness (W) increases from 0.822 to 0.851 when we add chemoattractant to the environment. Similarly, the most probable Holder exponent h 0 (q) shifts from 0.78 to 1.06 by adding the chemoattractant. This implies that in an environment rich in chemoattractant, the multi- fractal spectrum related to a bacterium motion gets a wider shape with higher values of h 0 (q), and this is true about all the density cases we analyzed except the highest density equal to 8x10 6 (bacteria/cm 3 ). In this case because of the very high density of bacteria in the environment, the volume exclusion effects and the receptor saturations there is no significant difference in h 0 (q) between the cases with and without chemoattractant and both of them are mono-fractal. In other words, adding chemoattractant to the environment impacts bacteria motion in a way that makes their motion more directional and less oscillatory. Figure 3.5c-d demonstrate the effect of obstacles in the environment for simulated E. coli. We compared the case of bacteria density 8x10 4 bacteria/cm 3 without chemoattractant and with obstacles (𝑊 = 1.004 and ℎ ` (𝑞)= 0.88) with the case of the same density and chemoattractant condition but without obstacles (𝑊 = 0.6918 and ℎ ` (𝑞)= 0.72). In an environment with obstacles the multi-fractal spectrum shifts toward lower values of h 0 (q) and narrower spectrum widths. In such an environment, the bacterium has the tendency to oscillate more than to follow a directional motion. 38 Figure 3.5 Multi-fractal analysis of motion for both simulated E. coli in BNSim and S. marcescens in vitro experiments (a) Multi- fractal spectrum of the simulated E. coli position in the direction of linear gradient of chemoattractant for different cases with and 39 without chemoattractant and different bacteria densities (case 2). (b) Generalized Hurst exponent of the simulated E. coli position in the direction of linear gradient of chemoattractant for different cases with and without chemoattractant and different bacteria density (case 2). (c) Multi-fractal spectrum of the simulated E. coli position in direction of linear gradient of chemoattractant with a density of 8x10 4 bacteria/cm 3 for different cases with and without chemoattractant and obstacles in the environment (case 1). (d) Generalized Hurst exponent of the simulated E. coli position in direction of linear gradient of chemoattractant with a density of 8x10 4 bacteria/cm 3 for different cases with and without chemoattractant and obstacles in the environment (case 1). (e) Multi-fractal spectrum of the experiments in vitro on S. marcescens position in the direction of linear gradient of chemoattractant for different cases with and without chemoattractant and different bacteria densities. (f) Generalized Hurst exponent of the experiments in vitro on S. marcescens position in the direction of linear gradient of chemoattractant for different cases with and without chemoattractant and different bacteria densities. (g) Multi-fractal spectrum of simulated S. marcescens position in the direction of linear gradient of chemoattractant for different cases with and without chemoattractant and different bacteria densities. (h) Generalized Hurst exponent of simulated S. marcescens position in the direction of linear gradient of chemoattractant for different cases with and without chemoattractant and different bacteria densities. 3.5 Investigation of Simulated S. marcescens Swimming Dynamics To investigate whether the mathematical characteristics observed for E. coli can also be retrieved from other bacteria, we modeled S. marcescens in BNSim (case 3) and extract the motion trajectories. By analyzing the results from this model we show simulated S. marcescens dynamics motion is highly nonlinear (Figure 3.6). Figure 3.6 Skewness plots for simulated s. marcescenns trajectories. Skewness plots for different bacteria density cases in simulation, with/out linear gradient of chemoattractants in y direction. In all the cases skewness is unequal to zero meaning bacterium displacement does not have normal distribution; therefore, bacterium motion is nonlinear. Figure 3.5g-h present the multi-fractal analysis for motion trajectories of simulated S. marcescens. The results show simulated S. marcescens bacterium exhibits multi-fractal dynamics. Increasing the bacteria population shifts the multi-fractal spectrum to the left. This implies that a single bacterium tends to show more mono-fractal type of behavior rather than multi-fractal type in a denser environment. The reasons of this behavior are the volume exclusion effects and chemical interactions between bacteria. Figure 3.5g shows that by adding chemoattractant to the environment multi-fractal spectrum shifts to the right for the Size of data [n] 0 2000 4000 6000 8000 10000 Skewness 0 100 200 300 400 500 600 700 800 10 4 (bacteria/cm 3 ), with Chemoattractants 10 4 (bacteria/cm 3 ), without Chemoattractants 10 5 (bacteria/cm 3 ), with Chemoattractants 10 5 (bacteria/cm 3 ), without Chemoattractants 10 6 (bacteria/cm 3 ), with Chemoattractants 10 6 (bacteria/cm 3 ), without Chemoattractants 10 7 (bacteria/cm 3 ), with Chemoattractants 10 7 (bacteria/cm 3 ), without Chemoattractants 40 same bacteria density. This implies that bacterium tends to move more directionally and oscillates less in the presence of chemoattractant. These results show that simulated S. marcescens and E. coli show similar mathematical characteristics. 3.6 Investigation of Real S. marcescens Swimming Dynamics Through In Vitro Experiments. We extract the S. marcescens motion trajectories from real experiments by single particle tracking method and analyze them. The results from skewness analysis show the same nonlinear behavior for S. marcescens motion (Figure 3.7). This is in agreement with the previously presented results for simulated E. coli and S. marcescens meaning that the motion of S. marcescens is highly nonlinear in the environment. Figure 3.7 Skewness plots for experimental S. marcescens trajectories. Skewness plots for different bacteria density cases in vivo environment, with/out linear gradient of chemoattractant in y direction. In all the cases skewness is unequal to zero meaning bacterium displacement does not have normal distribution; therefore, bacterium motion is nonlinear Figure 3.5e-f show the multi-fractal analysis for motion trajectories of S. marcescens from experiments in vitro. The results elucidate that S. marcescens has a multi-fractal motion in these environments. For the same population density by adding chemoattractant to the environment the spectrum moves towards right, which is a sign of multi-fractal behavior. We observe a similar but not identical multi-fractal behavior from real S. marcescens compared to simulated S. marcescens. The reason that the results do not match exactly is for the simulated S. marcescens we have control over most of the internal and environmental parameters. In contrast, real S. marcescens like other biological systems in nature show a higher degree of adaptation such as adjustment to environmental condition, self-organization and complexity in their behavior. Another reason for the difference between the simulation and experiment could be attributed to other forms of 50 100 150 200 250 300 350 400 0 50 100 150 200 250 300 Size of data [n] Skewness 10 6 (bacteria/cm 3 ), without chemoattractants 10 6 (bacteria/cm 3 ), with chemoattractants 10 7 (bacteria/cm 3 ), without chemoattractants 10 7 (bacteria/cm 3 ), with chemoattractants 10 8 (bacteria/cm 3 ), without chemoattractants 10 8 (bacteria/cm 3 ), with chemoattractants 41 interactions that happen in the biological world and are not captured in the current computational model. However, if someone discovers a new forms of interaction, they could be included in the computational model. This remains for the future work. 3.7 Summary Bacteria like E. coli and S. marcescens perform an efficient diffusive search (based on the status of their receptor signals) to find food in a complex environment [20]. In this chapter, we investigate the fractal characteristics of the bacterium motion while swimming in a three-dimensional environment. We first analyze the second and third order moment of single bacterium motion trajectory with the ergodicity and nonlinearity tests under different conditions including bacteria density and chemoattractant distribution. These help build the basis for higher order moment analysis of single bacterium motion through MFDFA method to capture the entire spatial and temporal complexity of single bacterium motion. Based on MFDFA analysis, we show bacteria swimming dynamics has multi-fractal characteristics. Moreover, we study how bacterial density and complexity of environment (i.e. obstacles and chemoattractant) influence the multi- fractal characteristics of dynamics of bacteria swimming. Consequently, our results contribute to the current state of knowledge by demonstrating that bacterial dynamics exhibits multi-fractal characteristics. The effect of chemoattractant is one of the main environmental conditions affecting bacteria multi- fractal characteristics. Multi-fractal analysis shows that in presence of an L-aspartate gradient in the environment, the multi-fractal spectrum gets wider and moves towards the right meaning that bacteria tend to move directionally and oscillate less compared to the case without a chemoattractant gradient being present in the environment. We also consider the effect of bacteria density as another important factor that affects the bacteria multi- fractal characteristics. Our results show that increasing the density of bacteria makes the multi-fractal spectrum of their motion shift toward the left and transit to a mono-fractal behavior. Moreover, further studies would be needed to relate these multi-fractal properties with the observed dependency between shear viscosity and density of the motile bacteria [177]. Previous research has shown that bacteria colonies display spatial fractal properties [178]. Ben-Jacob and his coworker studied the spatial fractal properties of Paenibacillus dedritiformis bacteria colonies under the effect of Septrin antibiotic [179]. They show in some patterns the members of the colony stay closer together and form large vortices, which can increase the colony’s ability to dilute the antibiotic with the lubricating fluid secreted by individual microbes; on the other hand, in some patterns their colony organizes into narrow and straight branches. This pattern happens when the food is scattered in the environment and will maximize bacteria contact with the limited nutrients in their environment. It is essential to emphasize 42 on the point that, in contrast to previous research, we investigate the multi-fractal characteristics of the bacterium motion while swimming in a three-dimensional environment in this study. This multi-fractal characteristic of the bacterium dynamics can get effected by continuous interactions between bacteria within the population. Moreover, this multi-fractal dynamic characterization of bacteria motion can allow us to quantify the effects of molecular (drugs) and gene-based therapies by investigating whether or not cells exhibit a directed movement as a result of the detection and transduction of external stimuli or continue to perform random wanderings. Interestingly, stem cells dynamics also exhibits a rich multi-fractal behavior, although stem cells do not use flagella or cilia for locomotion. This suggests that this mathematical framework can be possibly used for identifying universal mathematical characteristics for cell motion irrespective of the size, scale and nature of interactions [180]. A comprehensive understanding of individual and collective movement of cells can be important for developing mathematical models of cell growth that can abstract and account for biologically relevant features such as cell density and volume exclusion effects, mechanical forces acting on cells surface, environmental conditions (e.g., pH, temperature), chemical factors (e.g., chemo-attractants, chemo- repellants), as well as factors that affect the internal state (e.g., cell-cell communication [181], mitotic phase, cell health state, cell age). For instance, the above mentioned mathematical characterization can help decide the type of dynamical models (short-term memory vs. long-term memory), the level of detail at which individual cells must be described (how much of the internal state needs to be modeled), or the effect of various external stimuli on individual and collections of cells. Overall, this newly identified mathematical property can enable the microscopic description of the motion of individual cells through concepts and tools from non-equilibrium statistical physics (e.g., multi- fractal master equation). Integrating over the population density and the characterized effect of internal and external stimuli could enable us to construct macroscopic models and propose control strategies for the population of cells in complex environments to reach specific targets. An equally important research problem is to understand not only how the multi-fractal features allow for a compact multi-scale mathematical description, but also how it enables and simplifies the controllability of complex systems [182]. This is left for future work. 43 4 Mathematical Characterization of Bacteria Collective Motion 4.1 Introduction A group of agents or animals moving collectively is an example of a complex system. The collective characteristics between the agents are regulated by behavioral tendencies, as well as short-range and long- range interactions among them. The complex dynamics of such system evolves through different states (i.e., spatio-temporal arrangement of the agents inside the group, while moving collectively in the environment) [26], [27]. Despite significant research and progress in studying natural [28]–[36] and engineered [37]–[49] collective motion, the field is still trying to identify the dynamical states in a collective motion and predict the transition between them. Toward this end, we introduce the idea of using energy landscape to study the dynamics of a collective group (the content in this chapter have been published in [183]). We develop a new approach, which identifies and extracts the dynamical states of the spatial formation and structure for a collective group. Our mathematical framework enables the estimation of the free energy landscape [50]–[62] of the states of the group motion and also quantifies the transitions among them. In this approach, we are able to distinguish between stable and transition states according to their energy level and the length of time the group prefers to stay in each states. We noticed a collective group has a lower energy level at stable states compared to transition ones. This could be the reason that a group prefers to stay for a relatively longer time in stable states compared to transition states. Furthermore, the group’s structure may convert to transition states with higher energy level while reorganizing itself and evolving between two different stable states. The energy landscape helps to investigate the complexity of a collective motion. The remaining of this chapter is organized as follows: In section 4.2, we briefly overview the concept of energy landscape. In the next section, we present our algorithm to construct the energy landscape of a collective motion dynamics. Next, we construct the energy landscape of a simulated collective motion. In section 4.5, we quantify missing information of the group’s structure based on the energy landscape. Next, we quantify emergence and self-organization of a collective motion based on the missing information. In section 4.7, we estimate the complexity of the collective behavior. In section 4.8, we perform the complexity analysis based on energy landscape for three different natural collective behaviors. 44 4.2 Energy Landscape A complex system’s dynamics is influenced by interactions between the components of the systems. For example, consider the dynamics of a group of particles (e.g. atoms, molecules, proteins, etc.) as a complex system. In such systems, each element is coupled with the rest of the system based on the strength of the interactions between them. These interactions lead the dynamics of the system to evolve between different possible configurations defined as states of the system dynamics. For a complex system, some of the interactions between the elements of the system may be unknown. This makes it hard to predict and analyze such systems. Energy landscape can help to analyze such systems without a complete knowledge about the interactions between their components. Energy landscape provides a tool to identify and map all the possible states (i.e. configurations) of the system with their corresponding energy level (Gibbs free energy) [24], [25], [184]–[187]. From mathematical perspective, the energy landscape is a function representing an energy level to each possible state of the system (i.e. 𝑓: 𝑋 → 𝑅, where 𝑋 is a configuration for the system). If 𝑛 represents the number of degree of freedom of the system, we can consider 𝑋 = 𝑅 » . In this case, the plot representing energy landscape of the states is a hypersurface in 𝑅 »l a [188]. Figure 4.1 shows an example of a one degree of freedom energy landscape. The energy landscape represents all the possible states for a system. Hills and valleys in the energy landscape represents the states with local maxima and minima of the energy level 𝑓 respectively. Usually the states with the most populated entities in them corresponds to the regions of the map with lowest free energy (i.e. local minimums) [189]. These states can be considered as local equilibrium states of the system dynamics. Energy landscape provides a better understanding about the dynamics, structure and interactions between the elements of a complex system. It helps to identify the local equilibrium states of the system. Based on this landscape, the systems dynamics can be predictable relying on the initial condition of the system. In other words, knowing the initial condition of the system can make it possible to predict that the system will be absorbed to which local equilibrium state. Moreover, it can estimate the amount of required energy for a system to evolve between two different states on the landscape (Figure 4.1). Energy landscape can be defined in two level: (1) microscopic level, (2) macroscopic level. Microscopic level represents the dynamics of all the possible states of a system. On the other hand, in macroscopic level, all the possible states of the system will be partitioned into a limited number of macro- states. In this level we deal with the macro-states and their corresponding energy level. The free energy landscape concept was originally proposed to investigate potential energy in physical systems [190]. Free energy landscape concept developed more by Henry Eyring and Michael Polanyi during their research on the potential energy surface for the interactions between hydrogen atoms [191]. In 1980s, scientists started using energy landscape concept to study protein folding [192] and in 1990s the theory of 45 RNA free energy landscape have been developed [193]. Free energy landscape has been used as an approach to study the folding kinetics of RNA molecules [194]. From 1996, Energy landscape have been used to study the human brain through different activities [195]. Watanabe and his coworkers studied the dynamic of brain activity during bistable perception using free energy landscape [62]. In 2001, Adriaan de Lange developed a Gibbs free energy theory of human evolution. Baba and his coworker [54] developed an algorithm to extract the free energy landscape from a one-dimensional time series. Figure 4.1 Energy landscape of a complex system with one degree of freedom. 4.3 Algorithm to Construct the Energy Landscape of a Collective Motion We briefly overview our algorithm to identify and extract the dynamical states in a collective motion for a group of N agents moving in a three-dimensional space and build the energy landscape. Our framework is based on an approach presented by Akinori Baba and his coworkers [54], [196]. The steps are as follow: 1. We divide the time series containing the location of all the agents denoted by 𝑟(𝑡) to sub-intervals centered at time 𝑡 å with time window [𝑡 å −∆ 2,𝑡 å +∆ 2], where ∆ is the preferential time scale (Figure 4.2a). 2. We construct the probability density function of the location of all the agents in the group corresponding to each sub-interval (i.e. 𝑝 C ) and based on that we find cumulative distribution function (CDF) of the agents’ location in the space. We also estimate the CDF corresponding to the position for the entire group through the whole time in the same way. 3. Based on Kantorovich distance 𝑑 q we compare the CDF of sub-intervals with whole time series CDF and cluster the sub-intervals based on the similarities and closeness in the probability distribution function (equation (4.1)) [197]. Each cluster contains subintervals with similar dynamical configuration. 𝑑 q (𝑝 C ||𝑝 ¹ )= ( æ H® 𝑝 C 𝑟 d −𝑝 ¹ 𝑟 d ) 𝑑𝑟′ ® H® 𝑑𝑟 (4.1) 46 4. We consider each of the clusters as a spatio-temporal state for the group dynamics (Figure 4.2b). We calculate the escape time for each state, meaning the time between when the system enters and leaves each cluster. 5. We calculate the residential probability 𝑃 C of the ith state and transition probabilities 𝑃 C¹ from the ith state to the jth state (Figure 4.2c). 6. Based on these probabilities, we estimate the free energy landscape by quantifying the energy level in each state (F i ) from equation (4.2) and energy barrier for the group while evolving from state i to state j (F ij ) from equation (4.3) [196]. 𝐹 C = −𝑘 ] 𝑇 ln (𝑃 C ) (4.2) 𝐹 C¹ = −𝑘 ] 𝑇 ln ( ℎ 𝑘 ] 𝑇 𝑃 C¹ ) (4.3) In equation (4.2) and (4.3), symbol k B represents Boltzman constant. Symbols h and T are Plank constant and temperature respectively. Based on these energy levels we can estimate the free energy landscape of the group evolving between different states. Figure 4.2 Schematic description of the main steps for building the energy landscape for a group of N agents moving in a three- dimensional space. (a) First, we subdivide the trajectories of all agents in the group to equal sub-intervals centered at time 𝑡 å with a time window of [𝑡 å −∆ 2,𝑡 å +∆ 2], where ∆ is the predefined time scale. Next, we estimate the three-dimensional probability distribution function of the motion of the group for each sub-interval. (b) We use the Kantorovich metric to cluster these sub- interval time series based on their similarities in the probability distribution function. Each cluster of sub-intervals can be interpreted as a state for the collective motion. (c) In the last step, we estimate the transition probability matrix among the identified states of the collective motion. 4.4 Constructing the Energy Landscape of a Simulated Collective Motion To examine our framework, we first use a well-known agent-based model [37] which captures different behavior of a collection of interactive agents in a three dimensional space. This model is based on simplified local interactions between the individuals and is able to mimic the spatial dynamics of a group of animals such as bird flocks or fish school. By varying the degree of local interactions among the agents in this model, we can observe different types of behavior from the group [37]. The motion of each individual in the group is the outcome of local repulsion, alignment and attraction tendencies depending on the location 47 and orientation of the neighboring agents. The individuals tend to align themselves with the neighbors, while avoiding collision by keeping a minimum distance between them. Individuals avoid being isolated and keep the group to move as a single coherent entity by maintaining an attraction tendency between them. The details of the model are as follow. Consider N individuals in a group (i=1, 2, 3,…, N) with position vector 𝑝 C (t) and direction vector 𝑑 C (𝑡) at each time. The desired direction of each individual for the next time step based on local interaction between them is 𝑤 C (𝑡+𝜏) with 𝜏 representing the time step. This model considers three different zones around each individual (Figure 4.3). The first spherical zone called zone of repulsion and the individual is located at the center of it. If there are other neighbors in the zone of repulsion of an individual, the individual moves away from them to keep a minimum distance and prevent collision. The second spherical zone is the zone of orientation. If there is no other neighbor in the zone of repulsion of an individual, then the individual tries to align itself with other neighbors in its zone of orientation. The third spherical zone is the zone of attraction. The attraction between an individual and its neighbors in this zone results in the coherence of the group. Considering these three regions, there is a blind volume behind the individual in which the individual does not sense and respond to other neighbors in this zone. In this model, variables 𝑛 æ , 𝑛 ç and 𝑛 F represent correspondingly the number of neighbors in zone of repulsion, orientation and attraction of the agent. Variable 𝑤 æ (𝑡+𝜏) represents the desired direction of individual i with respect to repulsion from others in repulsion zone. 𝑤 æ 𝑡+𝜏 = − 𝑟 C¹ 𝑡 𝑟 C¹ 𝑡 » è ¹éC (4.4) 𝑟 C¹ 𝑡 = (𝑝 ¹ 𝑡 −𝑝 C 𝑡 ) 𝑝 ¹ 𝑡 −𝑝 C 𝑡 (4.5) When there are some neighbors in zone of repulsion (𝑛 æ ≠ 0), the individual i only reacts with respect to them. As a result, the desired direction 𝑤 C 𝑡+𝜏 = 𝑤 æ 𝑡+𝜏 can be quantified from equation (4.4) and equation (4.5). If there is no individual in the zone of repulsion, then the desired direction will be defined based on neighbors in zone of orientation and attraction (𝑤 C 𝑡+𝜏 = a e ×(𝑤 ç 𝑡+𝜏 +𝑤 F 𝑡+𝜏 )). 𝑤 ç 𝑡+𝜏 and 𝑤 F 𝑡+𝜏 can be quantified from equation Error! Reference source not found. and equation Error! Reference source not found.. 𝑤 ç 𝑡+𝜏 = 𝑑 ¹ 𝑡 𝑑 ¹ 𝑡 » ê ¹¼a (4.6) 48 𝑤 F 𝑡+𝜏 = 𝑟 C¹ 𝑡 𝑟 C¹ 𝑡 » m ¹éC (4.7) Considering the desired direction vector at each time step, if 𝑤 C 𝑡+𝜏 is less than maximum turning rate 𝜃, then 𝑑 C 𝑡+𝜏 = 𝑤 C 𝑡+𝜏 . On the other hand, if desired direction vector exceeds the maximum rate, then the individual rotates by angle of 𝜏×𝜃 towards the desired direction. Figure 4.3 Different zones of interaction around each individual in a group of agents moving in three-dimensional space in a model proposed by Couzin and his coworkers [37]: Zone of repulsion, zone of orientation and zone of attraction. There is a blind zone behind each individual, in which they do not sense and react towards other neighbors in that area. The dynamics of the group can change between four different common collective behaviors depending on the width of different zones around the individuals (Figure 4.4). These four collective motion behaviors identified by Couzin and coworkers [37] are: torus, swarm, dynamic parallel group and highly parallel group. The torus configuration emerges when the individuals rotate around an empty space. This happens when the zone of orientation is relatively small compared to zone of attraction. In this case individuals have a tiny zone of repulsion around them (Figure 4.4a). On the other hand, when individuals exhibit an attraction or repulsion behavior between themselves and there is no orientation behavior (consequently no parallel motion), the colony demonstrates a so-called swarm collective behavior (Figure 4.4b). The dynamic parallel group emerges by increasing the zone of orientation, which causes the individuals align with each other and makes the group more motile compared to the above-mentioned two cases (Figure 4.4c). The group shifts to a highly parallel group behavior when the zone of orientation is relatively bigger compared to the case of dynamic parallel group (Figure 4.4d). In this case, the individuals are in a highly aligned arrangement and the group is more motile compared to the previous cases [37]. 49 Figure 4.4 Various collective patterns of a simulated model of a group of agents moving in a three-dimensional space. (a) Torus: Individuals rotate around a center point within an empty space. (b) Swarm: Individuals show attraction and repulsion behavior between themselves and there is no orientation behavior consequently no parallel motion. (c) Dynamic parallel group: Individuals align with each other and make the group more motile compared to two previous cases. (d) Highly parallel group: Individuals are in a highly aligned arrangement and the group is more motile compared to previous cases. We analyze these four types of collective behavior separately using our free energy landscape framework and show that: (1) each behavior is a combination of various structural states and (2) the group is transitioning among these spatio-temporal states over time. We identify and extract these states as building blocks for each collective behavior type. Figure 4.5 summarizes the transition probabilities between them for each case. We use the estimated transition probabilities to compute the free energy landscape using equations (4.2) and (4.3). Comparison between different cases shows that swarm behavior evolves between more possible states and this confirms that there is lower level of arrangement in the group in swarm phase. As a result, the transition matrix and correspondingly its energy landscape for this case has more spikes and is less smooth compared to the other cases (Figure 4.5b). In contrast, for torus and dynamic parallel group behaviors there is more structural order due to alignment and parallel motion between individuals. Consequently, the collective dynamics of the group is characterized by fewer possible states when compared to swarm behavior and the transition matrix and energy landscape is smoother (Figure 4.5a,c). For highly parallel group all the individuals are completely aligned with each other and the structure of the group has the highest arrangement compared to the other behaviors, so the structure of the group evolves 50 between the lowest number of possible structural states and the transition probability matrix is less spiky (Figure 4.5d). Figure 4.5 Transition probabilities among the states identified in different collective behaviors of the simulated agent-based model. (a) Torus, the plot shows the transition probability between different states in this collective behavior. (b) Swarm, the group of agents has the highest number of states in this collective behavior and the landscape has more spikes and is less smooth compared to the other cases. (c) Dynamic parallel group, the transition probability looks similar to the torus phase, this similarity is due to preference of individuals to align their motion parallel to their neighbors. (d) Highly parallel group, the group has the lowest number of possible states in this phase and the landscape is less spiky due to high preference of individuals to move parallel with respect to each other. 4.5 Interpretation of Missing Information for a Collective Motion To further exploit this free energy landscape description, we develop an information theoretic framework for quantifying the degree of missing information of a collective group motion. In general, missing information can be defined as quantifiable structure or pattern in a system [44], [198], [199]. It can be used as a measure of internal order of a system and uncertainty. According to Shannon, missing information can be defined from equation (4.8). 𝐼 = − 𝑃 C C log𝑃 C (4.8) We define missing information for a collective motion as the level of missing communicated information between the agents due to their short-range and long-range interactions. This can be interpreted as the amount of information needed to specify the coupling between the agents and as a result the exact structural state of the group forming a specific spatial arrangement. 51 In our framework, we find the dominant states of a collective motion dynamics from energy landscape analysis. The energy landscape is based on the transition probability matrix between different states (Figure 4.2c). Based on the transition probability matrix we quantify: (1) missing information of each state (2) missing information of the entire motion of group combining all the possible states. To compute the missing information for each state, we find the probability transition matrix 𝑃 of the collective motion evolving from one state to the others. Then we find the missing information related to each state i (𝐼 C ) from this probability transition matrix using equation (4.9): 𝐼 C = − 𝑃 C,¹ ¹ log𝑃 C,¹ (4.9) Next, we find the missing information related to the entire group motion considering all the possible states. We define matrix Q containing the probability of all the possible cases of transitioning from one state to the other, independent of the initial state in the transition. The missing information for entire group motion (I) can be quantified by equation (4.10): 𝐼 = − 𝑄 C,¹ ¹ log𝑄 C,¹ C (4.10) It is important to emphasize the difference between matrix P from Q. In probability transition matrix P the sum of all the elements in each row is equal to one, while in probability matrix Q the sum of all the elements in the entire matrix is equal to one. The difference is due to the way of normalizing the probability matrix. 4.6 Quantifying Missing Information for a Simulated Collective Motion In what follows, we quantify the missing information of the group for different types of collective behavior related to the simulated model presented in section 4.4. The missing information has been defined as a way to quantify the number of possible different ways to arrange a specific system [63], [198]. Therefore, it estimates the level of disorder in a system. For a collective motion, we consider it as a level of missing communicated information between the agents due to their interactions with each other. This can be interpreted as the amount of information needed to specify the exact physical state of the group forming a specific spatial arrangement. This also means how much internal order or uncertainty a group formation has in a specific state. Figure 4.6a compares the missing information of the group for different collective behaviors with different number of individuals. In our analysis, we consider the same initial conditions for the individuals’ location and speed. The radii of individuals’ interaction zones generating similar collective behavior are identical irrespective of the population size of the group. Figure 4.6a shows the transition from a swarm phase to torus, dynamic parallel group and then highly parallel group; corresponding to this 52 transition, we observe that the missing information is decreasing due to more alignment among the neighbors and increasing internal order of the group structure. To further investigate the effect of local interactions on missing information related to the group structure, we increase the group density to 100 individuals, we fix the radius of zone of repulsion and zone of orientation, and we change the zone of attraction. We consider the same initial condition for the individuals’ location and speed in all the cases. We perform the same analysis and quantify the missing information of the collective group. Figure 4.6b (blue line) shows that by increasing the radius of zone of attraction while maintaining the same initial condition it results in decreasing the missing information. Increasing the zone of attraction causes the individuals to interact more with each other rather than being isolated and as a result the group tends to show more collective arrangement. This means there is more structural order and less missing information about the group formation. Similarly, fixing the radius of zone of repulsion and zone of attraction and increasing the radius of zone of orientation contributes to a reduction in the missing information (red line in Figure 4.6b). This means that the expansion of orientation zone makes the individuals to align with each other and move parallel, which implies an increasing degree of internal order within the group. Consequently, the missing information decreases. Figure 4.6 Quantifying the missing information of the entire simulated agent-based model for various interaction rules. (a) We quantify the missing information from the dynamics of a group of agents considering different interaction rules which causes various collective behavior in the group while considering the same initial condition for the agents position. This plot shows the transition from swarm phase to torus, dynamic parallel group and then highly parallel group and the fact that the missing information is decreasing due to an increase in the internal order of the group. (b) The quantified missing information extracted from the group dynamics when the population consists of 100 individuals for varying value of the radius of zone of orientation while radii of zone of repulsion and attraction is fixed (red line) and in other case varying value of the radius of zone of attraction while radii of zone of repulsion and orientation is fixed (blue line). 4.7 Interpretation of Emergence, Self-organization and Complexity for a Collective Motion Emergence. Based on the missing information corresponding to each identified state of a collective motion, we can define emergence for a collective motion. Emergence in a system generally refers to some 53 information or characteristics of a system that appear in some states while they are not present in other states of the system [198], [200], [201]. Emergence of a system can be quantified by equation (4.11). 𝐸 = 𝑎𝑚𝑜𝑢𝑡 𝑜𝑓 𝑚𝑖𝑠𝑠𝑖𝑛𝑔 𝑖𝑛𝑓𝑜𝑟𝑚𝑎𝑡𝑖𝑜𝑛 𝑙𝑜𝑠𝑡 𝑡ℎ𝑟𝑜𝑢𝑔ℎ 𝑡ℎ𝑒 𝑒𝑣𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝑜𝑓 𝑡ℎ𝑒 𝑠𝑦𝑠𝑡𝑒𝑚 𝑏𝑒𝑡𝑤𝑒𝑒𝑛 𝑡𝑤𝑜 𝑠𝑡𝑎𝑡𝑒𝑠 𝑎𝑚𝑜𝑢𝑛𝑡 𝑜𝑓 𝑚𝑖𝑠𝑠𝑖𝑛𝑔 𝑖𝑛𝑓𝑜𝑟𝑚𝑎𝑡𝑖𝑜𝑛 𝑖𝑛 𝑡ℎ𝑒 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 𝑠𝑡𝑎𝑡𝑒 𝑜𝑓 𝑡ℎ𝑒 𝑒𝑣𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝐸 = ð kñ Hð êòÄ ð kñ = 1− ð êòÄ ð kñ (4.11) For a collective motion, we define emergence as being proportional to the structural order gained by the system and quantified in statistical terms with respect to the ambiguity of the initial state. Therefore, higher emergence shows more interdependence in the group due to stronger interactions between the agents. This shows the transferred information between the agents results in the dependencies of their motion in the group. In our framework we find the relative emergence in each identified state with respect to first state with the highest probability and lowest level of energy in the landscape. Self-organization. Self-organization in a collective motion can be defined based on the missing information corresponding to each identified state. In general terms, we call a system self-organize when the internal dynamics of a system increases its organization in time [198]. equation (4.12) quantifies the self-organization of a system [198]. 𝑆 = 𝐼 C» − 𝐼 çón (4.12) For a collective motion, we consider self-organization as a measure for the transformed information between the agents into the internal order of the group structure. This could be considered as a measure for group intelligence. In our framework we find the level of increase in self-organization of the group going from any possible state to the first and most probable state to be able to compare the self-organization of different states in a collective motion. Complexity. In general, Complexity represents the balance between emergence (presents variety) and self-organization (presents order) of a system [198], [199], [202], [203]. In other words, it is a balance between ordered and chaotic dynamics. Complexity of a system can be quantified by equation (4.13) [198]. 𝐶 = 𝐸×𝑆 (4.13) We consider the complexity for a collective motion as a combination of interdependency between the agents and internal order in the group structure as a result of transferred information between the agents due to short-range and long-range interactions. In our framework, we quantify the relative complexity of each state in a collective motion with respect to the complexity of the first possible state from emergence and self- organization corresponding to that state. 54 4.8 Complexity Analysis for Three Natural Collective Motions Swimming Bacteria. Cellular collective groups can create complex patterns based on their simple behavioral rules. Different ways of communication and information transfer between them affects their individual behaviors and motion. To investigate and quantify the complexity of cellular groups, we study motion of S. marcescens bacteria with density of 10 8 (bacteria/cm 3 ) moving in vitro. We obtain the data for a group of S. marcescens moving in three-dimensional space from Edwards et al. [21] and Zhuang et al. [20]. We consider the motion of 9 bacteria (our analysis can be extended to a higher number if we are able to track and extract long time trajectories of all the individuals) in the experiments and consider it as a collective group for our analysis. Next, we identify the possible states for the spatial formation of the selected group by analyzing the time trajectories of individuals’ motion. We investigate two different cases with and without distribution of chemoattractant in the environment. Considering the case without chemoattractant distribution in the environment (Figure 4.7a), the analysis identifies four states for the formation of the group motion. When the transition probability of one state has higher level, it indicates that the amount of energy needed for the group of bacteria to create that formation is relatively lower. Based on our analysis, the first state has the lowest energy in the landscape compared to others. This state with the lower energy formation has the highest transition probability in the transition matrix, which means it is the most probable state for the group formation over time. We can order other states relatively to this state based on energy landscape from lowest energy to the highest. We categorize them into two groups of stable and transitioning states. When the group has higher probability to stay in the same state over time, we consider it as a stable state, which can be recognized as local equilibrium state for the group formation. These stable states have lower energy level compared to others. On the other hand, if it is not probable for the group to stay in the same state, we consider that as transient state. Because of the difference in the energy level, the group prefers to stay for a relatively longer time in stable states compared to transition states. Meanwhile, the group may shape the transition states while evolving between two different stable states. Figure 4.7a,c demonstrate that adding chemoattractant to the environment, decreases the number of possible states for the group dynamics and the free energy landscape is more smooth compared to the case without chemoattractant. The existence of chemoattractant in the environment contributes to the preferable alignment of the bacteria motion with each other and causes the group to get more organized with less oscillatory and scattered motions. 55 Figure 4.7 Transition probability matrix and complexity analysis for a collective group of bacteria. (a) Transition probabilities among the possible states for a collective group of 9 bacteria selected from a population density of 10 8 bacteria/cm 3 moving in an environment without chemoattractant gradient. (b) Complexity analysis for different states compared to the first identified state with lowest energy level for a group of 9 bacteria selected from a population density of 10 8 bacteria/cm 3 moving in an environment without chemoattractant gradient. This plot shows the level of change in missing information when the collective motion leaves each identified state to evolve to a new state. It also demonstrates the relative emergence and relative self-organization and relative complexity of the collective motion when evolving from any of the identified state to the first and most probable state. (c) Transition probabilities between the possible states for a group of 9 bacteria selected from a population density of 10 8 bacteria/cm 3 moving in an environment with chemoattractant gradient. (d) Complexity analysis for different states compared to the first identified state with lowest energy level for a group of 9 bacteria selected from a population density of 10 8 bacteria/cm 3 moving in an environment with chemoattractant gradient. Based on our analysis, we observe that the missing information of the first local equilibrium state is the lowest (Figure 4.7b,d). Accordingly, for the transient state the missing information is higher compared to the stable states. Figure 4.7 shows the level of change in missing information when the collective group leaves any of the identified states to evolve to a new state. Figure 4.8 makes the concept of level of change in missing information when the group motion leaving state i to any possible state more clear. In this figure, when the dynamic of the group leave state 1 and evolve to any other possible state, the level of missing information increases on average, irrespective of the terminal state. Therefore, we can conclude that the absolute missing information of state 1 compared to other states is lower. 56 Figure 4.8 Level of change in missing information through the evolution of the dynamic of collective group from state 1 to other states. We define emergence for a collective group as being proportional to the structural order gained by the system and quantified in statistical terms with respect to ambiguity of the initial state [200]. Based on this definition, changing the system dynamics from one state to another can be interpreted as a way to transform or propagate the information. Figure 4.8b,d show the relative emergence of all possible states with respect to the most probable state, which has the lowest amount of missing information. Based on our analysis, high emergence means more dependencies in the group as a result of stronger interactions. This indicates that the group has more interdependent components in stable states with lower level of missing information compared to transition states. On the other hand, low emergence shows the group has more independent components, which represents the group in transition states. We define a group to be self-organized when the internal dynamics of the group increases its organization over time. This could be a measure for collective group intelligence and we plan to investigate further in our future work. Our proposed free energy landscape analysis shows that the group is naturally attracted to be in stable states with lower energy compared to transitioning states with higher energy level. By considering these attractor stable states as organized states, the group will be self-organized through time. Figure 4.8b,d show the level of increase in self-organization of the group going from any possible state to the first stable state with lower energy level (i.e., each point shows the amount of increase in group self-organization evolving from the corresponding state to the first stable one). The figure shows the group gets more self-organized when it evolves from any possible state to the first and most probable one with the lowest level of energy. This shows that the bacteria group gets more self-organized over time. The final purpose of our analysis is studying the complexity of a group motion. To quantify the degree of complexity for a group with specific types and possibly unknown or impossible to detect agent-to-agent interactions, we compute a complexity metric as the product between emergence and self-organization. Figure 4.8b,d show the relative complexity of all the possible states with respect to the first and the most stable state. Each point in this plot shows how complexity changes by evolving from the corresponding 57 state (i.e., the states represented by that point) to the first stable one. This figure shows that the complexity metric exhibits an increasing tendency when the group evolves from transition states to stable ones. This shows that over time, the group tends to stay more in the stable states with higher complexity compared to other ones. Flying Pigeon. Next, we analyze two different types of flying pigeon groups: free flight and home flight. In free flight case, the pigeons are flying freely in the sky while in the home flight they are migrating from one region to another region. We obtain the data for group of Pigeons Flying in three-dimensional space from Nagy et al. [204]. Figure 4.9 show that for free flight we have more dominant states compared to the home flight. This demonstrates that when the group has a destination and its goal is to reach its destination rather than just flying freely in the sky, it oscillates between less number of dominant spatial state formations. Our analysis of the proposed information metric (see the results in Figure 4.9) demonstrates that the stable states have a lower degree of missing information and higher degree of emergence, self-organization and complexity compared to transition states. This means that over time, the group of pigeons, independent of their flight type, tends to have spatial formation / structure related to stable states which has lower energy, higher degree of complexity compared to the transition states. Figure 4.9 Transition probability matrix and complexity analysis for a collective group of pigeons. (a) Transition probabilities between the possible states for a group of 9 pigeons in free flight. (b) Complexity analysis for different states compared to the first identified state with lowest energy level for a group of 9 pigeons in free flight. (c) Transition probabilities between the possible 58 states for a group of 8 pigeons in home flight. (d) Complexity analysis for different states compared to the first identified state with lowest energy level for a group of 8 pigeons in home flight. Group of Ant. Insect’s societies can be considered as an example of complex systems. For instance, a group of ants exhibit emergent characteristic at a higher level compared to the sum of emergent corresponding to all individuals separately. This means the group reacts like a single coherent entity in different situations (e.g., presence of attack to different part of the group) [205]. Therefore, scientists consider a group of ants as a single super-organism [206]. The individual ants in a group tend to form spatial organized structure (i.e. spatio-temporal states) with respect to each other. Using our framework, we can identify these spatio-temporal states, build their energy landscape and quantify their complexity. Regarding this, we analyzed a group of eight ants with identical role inside their population with our algorithm. We obtain the data for ant trajectory from Queen Mary University database directory [207]. Figure 4.10 shows the transition probability matrix. In this figure, the high peak points correspond to the lower energy levels in the landscape, meaning that the transition of the group among these states consumes less energy. Figure 4.10 shows the missing information and complexity analysis. We can see the same pattern meaning that the stable states have lower missing information and higher emergence, self-organization and complexity compared to transition ones. Figure 4.10 Transition probability matrix and complexity analysis for a collective group of ant. (a) Transition probabilities between the possible states for a group of ants. (b) Complexity analysis for different states compared to the first identified state with lowest energy level for a group of ants. It has been shown that the spatial organization of ants with different duties in a group strongly depends on their role [205]. Consequently, depending on the roles they play in the group they may form different structural states. One of the potential applications of our framework in the future can be studying the performance of ants with the same role inside their colony under different environmental conditions (e.g., attack to different parts of the group, migrating to new nest). We can also quantify the information transfer among members of different species inside a colony and compare the dependency of communication 59 between them based on the role they are playing inside the group. These two potential applications of our framework remains for the future work due to lack of access to required data for these analyses. 4.9 Summary Animals moving in a group are influenced by their social context meaning they adjust their motion in response to interactions with their neighbors and environment [208]. They keep a minimum distance from each other to avoid collision. Meanwhile, they have a long range attraction to others, which keeps them united as a group and prevent their isolation from the rest of the group. At the same time, they tend to align the direction of their motion with the ones near them to move in a synchronous fashion. These interactions between agents in the group are due to their sensory systems including vision, smell detection / chemical processing and sound. These multi-modal heterogeneous interactions among the agents cause the motion of the group to evolve through various spatio-temporal structures while moving as a synchronized and coherent entity without having a centralized controller. Such synchrony and structural patterns in the group helps the individuals to amplify their sensitivity and reactions / agility to the environmental conditions while they have limited individual sensing and processing capabilities. This proves critical for their survival [208], [209]. As an example, when a predator attacks the group, a small portion of the group senses the attack prior to the rest of the group. The efficiency of achieving a high degree of collective behavior helps to adapt faster to perturbation and decreases the reaction time of the whole group to the dangerous situations. In this case, they transform to a specific structural pattern to align more strongly with each other helping them to escape faster from the threat. As another example, such synchrony between them helps the group to identify the resources in the environment more efficiently. Hence, the spatial structuring within groups has important and evolutionary consequences [37], [210]. From this perspective, it is crucial to be able to study the whole group considering their structure and to develop a mathematical framework for identifying and quantifying the information flow within the group. Our mathematical framework helps to analyze various types of collective behavior exhibited by a group and identify / extract the possible spatio-temporal states that correspond to the highly interdependent group dynamics. Estimating the transition probability matrix between different states helps us to construct the energy landscape of the collective group evolving among these possible states over time. Relaying on this probabilistic characterization of the interdependency structure among various states, we quantify the missing information corresponding to the structural formation of a particular collective. We show that when local interactions among individuals increase in strength, the individuals tend to align more with their 60 neighbors and as a result the group gains more internal order. Therefore, the group structure does not change too much through time and as a result the number of possible states decreases and the missing information of the group structure decreases as well. We believe that this will help us understand how group of moving agents overcome the information bottleneck and plan to design new real experiments [211]. We also quantify the missing information, emergence, self-organization and complexity of the group corresponding to each of its possible structural states. We show that over time the group tends to stay in stable states with lower level of energy; this corresponds to higher degree of self-organization and complexity compared to other possible states. Our analysis demonstrates that the complexity of the group formation increases over time, which could be attributed to the fact that the interactions are evolving or adapting to external cues. Our mathematical framework can help us understand the evolution of behavior of various complex systems, from human microbiome to road traffic and potentially also economic and social networks. An important applicability domain of the proposed framework is represented by the need for a robust mathematical formalism for quantifying the efficiency, adaptivity, robustness and agility of a group of artificial learning cells and comparing how two artificial groups with different heterogeneous interactions and learning capabilities can perform on different environments with various degrees of uncertainty. Our framework could also serve as an initial step towards a solution to one of the main challenges in collective motion optimization and control. Communication between agents enables a decentralized control strategy for collective motion optimization. This causes the group of agents to self-organize and creates spatio-temporal patterns and ordered structures while following a good path at a specific time for their motion. This optimization observed in the group motion is a sign of intelligent behavior. According to Gerardo Beni [209] an intelligent group can be considered as a large parallel computational system, which performs computation and motion in parallel. Computation and simulation time for an agent based model, which predicts the group performance from its initial state scales with the number of agents. If the number of agents is large enough then the computation time increases exponentially and also the possible outcome after a certain finite number of steps of evolution of the group is a NP-complete problem [38], [47], [212], [213]. Therefore, control of such group with decentralized controllers is still a fundamental challenge, because there is often no obvious relation between the individual’s behavior and the final behavior of the whole group [39]. Our algorithmic strategy can be integrated into an engineering framework to be used to set the parameters that governs the dynamic of one agent and its corresponding interactions with other agents to achieve a certain degree of emergence, self-organization and complexity (Figure 4.11). For instance, one could combine an agent based model [28], [37] (where the characteristics and behavior of each individual are driven by several knobs) and macroscopic analytical model [47], [214] to describe the 61 collaboration in a group of reactive robots. In the next step, use our framework and quantify the emergence, self-organization and complexity as a function of the control knobs and size of the group. This combined analysis can enable the identification of the critical design considerations. Consequently, by playing with the local interactions between the agents, we can regulate the system to evolve towards desired states and control the corresponding free energy landscape. Controlling over the energy landscape implies following a few rules of interaction that contribute to a particular set of states with desired degree of self-organization and complexity. This remains for future work. Using our framework can help to build, characterize and optimize a group of robots with much simpler components with a decentralized control characteristic and replace the complex centralized control systems to perform the same task. By comparison, simplicity of the agents and decentralized control characteristics of the group make it possible for the group to adapt dynamically better to the environment and recover from different disturbances in the environment [209]. Therefore, the collective group systems can be more reliable to survive through disturbance compared to centralized control systems. Figure 4.11 The proposed algorithm can be integrated into an engineered framework which could help to set the parameters representing the dynamics of the agents and their interactions to achieve a certain degree of emergence, self-organization and complexity. 62 5 Mathematical Analysis of Bacteria Swarms for Fighting Infections: A Case Study on Ebola Treatment Approaches 5.1 Introduction Ebola outbreak in 2014 is one of the recent largest epidemics in history with devastating consequences on human society. Ebola has been found for the first time in central Africa in 1976. The virus survived inside a type of bat and human got infected for the first time by contact with infected bat. The virus spread through contaminated body fluids from one individual to another. Some of the symptoms of the infection are sever fever, diarrhea, vomiting, and 70% of the infected people die seven days after the symptoms appear [215]. When the virus enters human body, immune system response acts immediately by neutralizing the virus and trying to clear the body from the infection [216]. At first, the virus will be examined by the macrophages to be identified for human body immune system. This will help immune system to initiate specific defense mechanism against the infection. Then the macrophages present their findings to the CD4 positive T lymphocytes (CD4+T cells). These cells are known as helper T cells, which serve as the command center for the immune system [217]. The “CD4” represent a protein marker on the surface of the T cell. The “T” relates to thymus, which is an organ maturating the cell after being manufactured and migrated from bone marrow [217]. Depending on the information provided for T cells, if the immune system needs to respond, the helper T cells reproduce to increase command forces. Moreover, immune system also activates a second type of T cell, the CD8 positive T lymphocytes (CD8+ T cells). These T cells, also known as killer T cells, destroy infected cells with pathogens [217], [218]. Mathematical model of the dynamical interactions between the Ebola virions and the human immune system at cellular level can provide a comprehensive understanding of how the immune system functions during an Ebola infection. This can provide a toll to study the efficiency of different treatment strategies to cure the infection. Related to this, Herz presented a mathematical model for viral dynamics [219]. Later, rational vaccine design model have been proposed for Ebola virus based on Herz model [220]. These models still remain imprecise, because they lack the self-proliferation of T cells in presence of virus. They also do not consider the growth of virus from other type of infected cells in the environment. Considering these facts, in this chapter, we develop a new mathematical model of the dynamical interactions between the Ebola virions and the human immune system at cellular level. This model can 63 provide a comprehensive understanding of how the immune system functions during an Ebola infection. Next, we present mathematical models for three potential treatments for Ebola infection: chemotherapy, gene therapy, and bacteria therapy. Chemotherapy and gene therapy are two treatment approaches for other similar viral infection (i.e., HIV), which has almost identical dynamic behavior as Ebola virus while interacting with immune system [220], [221]. During chemotherapy, a specific dose of drug will be delivered to infected human body periodically, which reduces viral replication [217], [222]–[224]. On the other hand, in gene therapy, engineered genes are delivered ex vivo, and potentially ultimately in vivo, into a patient’s cells [225], [226]. They act within these protected cells to disrupt the viral life cycle. Unlike two previous approaches, in bacteria therapy, engineered bacteria will be injected to the body, they will colonize and secrete a chemical as a result of quorum sensing [227]. This chemical phosphorylates and deactivates Cathepsins enzymes in immune system, which is required for Ebola attachment and entry into human T cells [228]–[230]. As a result, it will prevent infection of T cells. The models presented in this chapter offers a precise paradigm to study the efficiency of different treatment strategies to cure the infection. This helps us to overcome drug design limitation and minimize the trial approach needed to design the drug while allowing to take into account the personal characteristics of a patient through precise parameter tuning. Moreover, we also quantify the efficiency of chemotherapy, gene therapy, and bacteria therapy approaches to cure Ebola infection. Overall, our new results demonstrate that bacteria therapy and gene therapy can be two applicable strategies to overcome the progress of Ebola infection inside host body. In previous chapters of this thesis, we investigated mathematical modeling of bacteria motion both at cell level and swarm level. Our motivation to do this study was using the bacteria for therapeutic purposes. Based on this, we propose using bacteria to treat infectious disease like Ebola in this chapter. The remaining of the chapter is organized as follows: Section 5.2 provides an overview of the Ebola virus dynamics while infecting the human cells. Section 5.3 presents our mathematical modeling for a network of healthy T cells and infected 𝑇 C cells and virus entities. We also present our model for three different treatment approaches. Section 5.7 describes the strategy for estimating the model parameters to match the real dynamic behavior of the virus during human body infection. Section 5.8 summarizes the analysis of our model and presents the fluctuations of different species during the treatment. 5.2 Ebola Virus Dynamics While Infecting the Human Cells Ebola does not have the ability to reproduce independently. Therefore, it must find a suitable host to provide the virus a comfortable environment to reproduce in a sequential five-step process: connection to the host cell, penetration, replication of the virus inside the host cell, assembly of the new viruses, and burst of host cell to release the viruses (Figure 5.1) [220]. Ebola virus targets CD4+T cells, which play a key role in the 64 immune system of human body. The virus carries a copy of its RNA. When the body is infected, the virus attaches itself to the T cell and injects its content, including the copy of its RNA, into the host T cell. Then the RNA will be transcribed into DNA by an enzyme called reverse transcriptase, which carried and injected to the host cell by virus. At this step the DNA of the virus will be duplicated inside the host cell and incorporated into the host genome. Next, new viruses sprout from the surface of the host cell by sparing or bursting the host cell [225]. Figure 5.1 Diagram of Ebola infection. Cell level interaction of Ebola virus with the human immune system T cells. Virus inters human T cells and infects them. Bursting infected cells release copies of virus. CD8+T cells are another class of T cells in immune system, which respond to infection. These cells cannot get infected with the virus, but they are able to destroy infected CD4+T cells. Therefore these cells in the immune system are responsible to kill the virus replicated inside the infected T cells [217] . In our modeling framework, the T cells include the two types of CD4+T and CD8+T cells. 5.3 Mathematical Model Overview In what follows, we develop a mean-field model describing the interaction of the Ebola virus with human immune system cells. Our model helps to monitor the fluctuations within the virus and immune system cells populations while they interact during infection. This model encapsulates all events in viral reproduction using a network of three entities: T(t) is the population of healthy target CD4+T cells; 𝑇 C 𝑡 is the population of productively infected T cells; and V(t) is the population of free virus. We assume that all the equations are scaled appropriately meaning that all the cell populations in our model occupy the same environment size. The model is as follows: Equation (5.1) describes changes in the population of healthy T cells. The first term in this equation shows the thymus-input source for new T cells. The next term represents the natural 65 death of T cells because of their average finite life span equal to a ô õ . Third term shows the self-proliferation of T cells in presence of the virus. The variable r represents the maximal proliferation rate, and h stands for half saturation constant of the proliferation process. The last term shows infection of T cells by virus with constant rate of k (when there is no treatment d(t)=1). Equation (5.2) shows the fluctuation in the population of infected 𝑇 C cells. The first term represents a gain for 𝑇 C from infection of T cells by the virus. The next term shows the loss of infected cells by their natural finite life span with rate of 𝜇 º k. The last term shows the loss of infected cells because of being stimulated to proliferate. During the proliferation the infected cells burst because of the large viral load. Equation (5.3) represents free virus population dynamics. First term is the source of virus population released by the burst of infected cells. N is the average number of virions released from every single infected cell. The second term shows loss of virus by the specific immune response in which the CD8 + T cells kill virus with rate of k T . The last term shows growth of virus from other type of infected cells in the environment with growth rate of g V and the half saturation constant b. 𝑑𝑇 𝑡 𝑑𝑡 = 𝑠 𝑡 − 𝜇 º 𝑇 𝑡 + 𝑟 𝑇 𝑡 ℎ+𝑇 𝑡 −𝑑 𝑡 . 𝑘 𝑇 𝑡 𝑉 𝑡 (5.1) 𝑑𝑇 C (𝑡) 𝑑𝑡 = 𝑑 𝑡 . 𝑘 𝑇 𝑡 𝑉 𝑡 −𝜇 º k𝑇 C 𝑡 − 𝑟 𝑇 C 𝑡 ℎ+𝑇 C 𝑡 (5.2) 𝑑𝑉(𝑡) 𝑑𝑡 =𝑁𝑟 𝑇 C 𝑡 −𝑘 º 𝑉 𝑡 + 𝑔 ö 𝑉 𝑡 𝑏+𝑉 𝑡 (5.3) Figure 5.3 in Section 5.7 shows the fluctuation of all the species in the Ebola infection network based on analyzing our model, which reproduces the reality. Based on our model the population of healthy T cells drops below 10% of its initial population, which represents a dying patient. In the next step we study how different treatment strategies affect the dynamics of this viral network. This will help us to optimize the process of drug design for Ebola infection. 5.4 Virus Infection Model Resistance to Chemotherapy One possible way of virus treatment is chemotherapy. In this approach, drugs work as inhibitors of reverse transcriptase. The role of such drug is to restrict the transcription of RNA to DNA. Therefore, it will stop the infection process and control the viral spread. Unfortunately, such drugs are not able to cure the infection. Instead they are able to control the exponential progression of the virus [217]. To include chemotherapy in our model, we multiplied parameter rate k by a function d(t) which is “on” when there is treatment and “off” when the system is outside treatment period [217]. 66 𝑑 𝑡 = 1 𝑜𝑢𝑡𝑠𝑖𝑑𝑒 𝑡ℎ𝑒 𝑡𝑟𝑒𝑎𝑡𝑚𝑒𝑛𝑡 𝑝𝑒𝑟𝑖𝑜𝑑,"off" 𝑝ℎ𝑎𝑠𝑒 𝑒(𝑡) 𝑑𝑢𝑟𝑖𝑛𝑔 𝑡ℎ𝑒 𝑡𝑟𝑒𝑎𝑡𝑚𝑒𝑛𝑡 𝑝𝑒𝑟𝑖𝑜𝑑,"on" 𝑝ℎ𝑎𝑠𝑒 In the “on” phase of treatment the viral infectivity is decreased which we model by defining e(t) as a treatment function with the following boundaries: 0 < e(t) < 1 presenting the percent effectiveness of the drug. Chemotherapy is a dose dependent treatment. So, e(t) can be considered as a factor representing the dose and also efficiency of the drug. Smaller e(t) represents a drug with higher efficiency. Outside the treatment period (i.e. “off” phase) we defined function d(t) = 1. This model helps us to study the effect of a variety of drugs with different efficiencies. We can also find an optimized treatment period for the chemotherapy. Such mathematical model can help us to reduce the trial process of designing the drug. 5.5 Virus Infection Model Resistance to Gene Therapy In order to account for a gene therapy treatment in our proposed model, we added two more cell types to our network besides T, 𝑇 C , and V entities: P and 𝑃 C . We consider two types of susceptible CD4+T cells: Protected P cells are type of transduced cells infused to the patient body for treatment, and natural T cells that were not manipulated and we introduced it in the last model. The 𝑃 C represents productively infected P cells. Figure 5.2 shows the network diagram characterizing the interactions between the Ebola virus and the immune system cells during gene therapy. Equation (5.6) and (5.7) track the changes in the population of P and 𝑃 C cells in our model. Figure 5.2 Diagram of gene therapy for Ebola infection. Cell level dynamics interactions of Ebola virus and the human immune system T cells during gene therapy. Ebola virus can enter both T cells and P cells and infect them. 67 The first term in (5.6) shows self-proliferation of P cells in presence of virus. The next term represents the natural death of cells because of their average finite life span with rate of 𝜇 º . The last term shows infection of P cells by virus with constant rate of 𝜀k. In (5.7), the first term represents a gain for 𝑃 C from infection of P cells by virus. The next term shows the loss of infected cells by their natural finite life span with rate of 𝜇 º k. The last term shows the loss of infected cells because of being stimulated to proliferate. 𝑑𝑇 𝑡 𝑑𝑡 = 𝑠 𝑡 − 𝜇 º 𝑇 𝑡 + 𝑟 𝑇 𝑡 ℎ+𝑇 𝑡 +𝑃 𝑡 −𝑘 𝑇 𝑡 𝑉 𝑡 (5.4) 𝑑𝑇 C (𝑡) 𝑑𝑡 = 𝑘 𝑇 𝑡 𝑉 𝑡 −𝜇 º k𝑇 C 𝑡 − 𝑟 𝑇 C 𝑡 ℎ+𝑇 C 𝑡 +𝑃 C 𝑡 (5.5) 𝑑𝑃 𝑡 𝑑𝑡 = 𝑟 𝑃 𝑡 ℎ+𝑇 𝑡 +𝑃(𝑡) − 𝜇 º 𝑃 𝑡 −𝜀 𝑘 𝑃 𝑡 𝑉 𝑡 (5.6) 𝑑𝑃 C (𝑡) 𝑑𝑡 = 𝜀 𝑘 𝑃 𝑡 𝑉 𝑡 −𝜇 º k𝑃 C 𝑡 − 𝑟 𝑃 C 𝑡 ℎ+𝑇 C 𝑡 +𝑃 C 𝑡 (5.7) 𝑑𝑉(𝑡) 𝑑𝑡 =𝑁𝑟 (𝑇 C 𝑡 +𝑃 C 𝑡 )−𝑘 º 𝑉 𝑡 + 𝑔 ö 𝑉 𝑡 𝑏+𝑉 𝑡 (5.8) 5.6 Virus Infection Model Resistance to Bacteria Therapy We propose bacteria therapy for Ebola treatment as a potential approach. Ebola virus needs to attach itself to human immune system cells through protein bonds to be able to transfer the copy of its RNA to the cell and infect it [202], [228], [231]–[233]. Cathepsins B and L are two required proteases in human immune system cells for Ebola attachments and entry to the cell [234]. In this therapeutic strategy, the engineered bacteria will be injected to the patient body. The bacteria population will grow inside the body. Quorum sensing between bacteria in the swarm causes releasing a chemical, which is needed for phosphorylation and deactivation of Cathepsin protein. This will prevent Ebola virus from attaching and entering immune system cells. This approach can help to prevent the progress of infection inside the host body. To consider bacteria therapy treatment in our proposed model, we added two more species to our network besides T, 𝑇 C , and V entities: M and B. M is the population of modified protected cells from attaching Ebola virus which cannot be infected by the virus and B in the population of engineered bacteria. Equation (5.11) and (5.12) track the changes in the population of M and B cells in our model. The last term in (5.9) represents the transform of T cells to M cells by deactivation of Cathepsin with rate of 𝐾 ] . The first term in (5.11) shows self-proliferation of M cells in presence of virus. The next term represents the natural death of cells because of their average finite life span with rate of 𝜇 º . The last term in represents the 68 transform of T cells to M cells by deactivation of Cathepsin. First term in (5.12) shows the natural growth rate of bacteria inside body and the second term represents the death of these bacteria with rate of 𝜇 ] . 𝑑𝑇 𝑡 𝑑𝑡 = 𝑠 𝑡 − 𝜇 º 𝑇 𝑡 + 𝑟 𝑇 𝑡 ℎ+𝑇 𝑡 +𝑀 𝑡 − 𝑘 𝑇 𝑡 𝑉 𝑡 − 𝑘 ] 𝑇 𝑡 𝐵(𝑡) (5.9) 𝑑𝑇 C (𝑡) 𝑑𝑡 = 𝑘 𝑇 𝑡 𝑉 𝑡 −𝜇 º k𝑇 C 𝑡 − 𝑟 𝑇 C 𝑡 ℎ+𝑇 C 𝑡 (5.10) 𝑑𝑀 𝑡 𝑑𝑡 = 𝑟 𝑀 𝑡 ℎ+𝑇 𝑡 +𝑀(𝑡) − 𝜇 º 𝑀 𝑡 +𝑘 ] 𝑇 𝑡 𝐵 𝑡 (5.11) 𝑑𝐵(𝑡) 𝑑𝑡 = 𝑒 𝐵(𝑡) ℎ+𝐵(𝑡) −𝜇 ] 𝐵 𝑡 (5.12) 𝑑𝑉 𝑡 𝑑𝑡 =𝑁𝑟 𝑇 C 𝑡 −𝑘 º 𝑉 𝑡 + 𝑔 ö 𝑉 𝑡 𝑏+𝑉 𝑡 (5.13) 5.7 System Identification When Ebola virus enters human body, even at a low proportion of free virus particles to uninfected cells equal to 0.01, virions are capable of destroying 80-90% of targeted immune system T cells in seven days. We estimated the parameters of our model without treatment (we consider d(t)=1 for the case without treatment in (5.1)-(5.3).) to match this realistic assumption. Table 5.1 presents the parameters estimated for the fitted model. Figure 5.3 shows the population dynamic of different species in our model based on validated model. Figure 5.3a exhibits that the population of targeted T cells reaches less than 10% of its initial value which indicates the collapse of immune system in human body and as a consequence death of the infected patient. In the next step, we study our validated model under three different mentioned treatment strategies. 69 Table 5.1 Parameters and variables of Ebola model Symbol Variable/Parameter Value T Healthy target CD4+T cell population 𝑇 ` = [c/𝜇𝐿] 𝑇 C Infected CD4+T cell population 𝑇 ` C = 0 [c/𝜇𝐿] V Ebola free virus population 𝑉 ` = 10 [c/𝜇𝐿] P Healthy protected cell population 𝑃 ` = 10 [c/𝜇𝐿] 𝑃 C M B Infected protected cell population Modified protected healthy cell population Engineered bacteria population 𝑃 ` C = 0 [c/𝜇𝐿] 𝑀 ` = 0 [c/𝜇𝐿] 𝐵 ` = 10 [c/𝜇𝐿] s(t) Thymus-input source for new T cells 0.2 [c/𝜇𝐿×𝑑] 𝜇 º Death rate of uninfected T cell population 0.015 [1/𝑑] r Maximal proliferation of T cell population 20 [c/𝜇𝐿×𝑑] h Half saturation constant of the proliferation process 90 [c/𝜇𝐿] k Rate T cells become infected by free virus 0.00008 [𝜇𝐿/𝑐×𝑑] 𝜇 º k Death rate of infected T cell population 0.5 [1/𝑑] 𝜀 Inhibition factor 0.01 N Number of average free virus released from bursting an infected cell 12 𝑔 o Growth rate of virus from other sources than T cells 2 [c/𝜇𝐿×𝑑] b h e 𝜇 ] 𝑘 ] Half saturation constant of the external virus source Half saturation constant of the bacteria source Maximal proliferation of B cell population Death rate of bacteria population Rate T cells become protected 10 [c/𝜇𝐿] 80 [c/𝜇𝐿] 3×12 [c/𝜇𝐿×𝑑] 0.02 [1/𝑑] 0.008 [𝜇𝐿/𝑐×𝑑] c = cell, d = day, L = liter, 𝜇 = micro 70 Figure 5.3 Fluctuation of different species in the cell level network of Ebola infection. a) Population dynamics of healthy and infected T cells during Ebola infection. Healthy T cells population reaches less than 10% of its initial population after seven days. b) Growth of virus in the first ten days of infection. 5.8 Simulation Results and Analysis We analyzed the network of healthy T cells, infected cells, and virus during chemotherapy for different treatment periods ∆ and different drug efficiency e(t). Figure 5.4 represents the dynamics of T cell and virus population during the treatment for different cases. Our analysis covers a range of different treatment period ∆ with the same drug efficiency e(t) to study the effect of ∆ on chemotherapy efficiency. We also considered different cases with the same ∆ and different drug efficiency e(t). Figure 5.4a shows that, independent of treatment period ∆, nearly after 10 days the population of healthy T cells in the network reaches less that 20% of its initial value resulting in the collapse of immune system and death of the patient. This motivated us to consider the ideal case of constant rate of most efficient drug treatment, meaning e(t) is the smallest possible values and stays constant in this case. Even in this ideal case of chemotherapy the population of healthy T cells reaches 20% of its initial value meaning that the treatment was not efficient to prevent progress of the infection. Comparing different drugs with different efficiency e(t) with the same treatment period ∆ exhibits that the more efficient drug meaning smaller value of e(t) result in less decrease in the population of healthy T cells (Figure 5.4a). Figure 5.4b demonstrates the fluctuation in the population of Ebola virus during the treatment. For all the cases after a couple of days the population of virus reaches nearly 10 4 times of its initial population which shows the fast progression of infection inside body. Our results show that chemotherapy may not be an efficient treatment strategy to control and cure the Ebola infection. Hence, we need to consider more modern treatment techniques like gene therapy. 71 Figure 5.4 Population fluctuation of cell level network of immune system T cells and Ebola virus during chemotherapy. Symbol ∆ represents the treatment period and P is the efficiency of the drug during treatment. For different cases of drug efficiency with various treatment period the T cell population reaches 20% or less than its initial population in the first 30 days of Ebola infection. In the next step, we analyzed the network of T cells and virus considering protected cells in the network due to gene therapy treatment. Figure 5.5 shows the changes in population of both protected and non- protected T cells and infected cells and virus during gene therapy. This figure presents that during the treatment the population of non-protected T cells decreases. In contrast, as a result of self-proliferation in the existence of virus, the population of protected T cell increases in the network. Figure 5.5 Fluctuation of different species in the cell level network of Ebola infection during gene therapy. a) Population dynamics of healthy and infected T cells and also healthy protected P cells during Ebola infection. b) Fluctuation of infected protected P cells c) Growth of virus under gene therapy in the first forty days of infection. 72 To compare the efficiency of gene therapy based on protected cell population, we considered different ratio of unprotected T cell and protected ones in the beginning of infection. Figure 5.6 shows that decreasing the ration of the target T cell to protected cell population results in less loss in T cells population and counted as more efficient treatment. Figure 5.6b demonstrates that increasing the number of initial protected cells results in limiting the viral infection inside the body. Figure 5.6 Population fluctuations of different species in the network of immune system T cells, Protected P cells and Ebola virus during gene therapy. T 0 is the initial population of unprotected CD4+T cells, and P 0 is the initial population of protected cells. T cells include both protected and unprotected cells population. Based on these results, decreasing the ratio of initial unprotected to protected T cells results in less reduction in T cell population, which shows higher efficiency of gene therapy. In the last step, we analyzed the network of T cells and virus under bacteria therapy treatment started in the second day of infection. Figure 5.7 shows the changes in population of both infected and non-infected cells and virus during bacteria therapy. This figure presents the population of healthy cells decrease by %20 percent of its initial value in the first 10 days and then the engineered bacteria swarm affects the network and help immune system to recover the loss of healthy T cell. In this network, the engineered bacteria swarm restricts the growth of virus and the healthy immune system T cells are able to kill the entire virus population after 20 days. 73 Figure 5.7 Fluctuation of different species in the cell level network of Ebola infection during bacteria therapy. a) Population dynamics of healthy cells during Ebola infection. b) Fluctuation of infected T cells c) Growth of virus under bacteria therapy in the first forty days of infection. Based on our analysis, both gene therapy and bacteria therapy approaches are promising strategy to cure the Ebola infection. However, bacteria therapy is the most efficient approach recommended for Ebola infection treatment by limiting the loss of healthy cells in immune system (Figure 5.8). Figure 5.8 Comparison of three different therapeutic strategies. Bacteria therapy has less damage on uninfected cell population and it will cure the infection in the shortest amount of time. 74 6 Conclusion and Future Research Directions 6.1 Major Contributions of This Dissertation One of the main challenges in modern medicine is the detection of silent progression and migration of various disease inside the human body [106]. To overcome these challenges, we can rely on bacteria motility and engineer bacteria swarms capable of detection of different disease in early stages and targeted drug delivery inside the human body [235]. With the goal of using bacteria for therapeutic purposes, one of the main steps is to be able to predict their motion in the environment and be able to model it with mathematical tools. It is important to learn the nature of single bacterium motion and how it gets affected from environmental condition and interactions with the neighbors. Toward this goal, in this thesis, we focused on studying the motion of bacteria in two different scales of cell-level and population-level. Subsequently, we summarize our main contributions: • Traditional modeling of bacteria populations relies on simplifying assumptions such as the memoryless dynamics of their motion and ignores most of the physical and chemical interactions among them. To improve this state of affairs, in chapter 3 of this thesis, we study the dynamics of populations of bacteria both in vitro and in silico, while accounting for realistic conditions like volume exclusion, chemical interactions, obstacles, and distribution of chemoattractant in the environment. We demonstrate that, more than nonergodicity, the bacterial swarm dynamics is nonlinear and multi-fractal in nature. Finally, we demonstrate that bacteria dynamic behavior strongly depends on bacteria density and chemoattractant concentration. These mathematical characteristics represent the canonical set of features that need to be considered for the accurate modeling of bacteria motion based on the nature of their swarm dynamics. • Biological systems are frequently categorized as complex systems due to their capabilities of generating spatio-temporal structures from apparent random decisions. In spite of research on analyzing biological systems, we lack a quantifiable framework for measuring their complexity. To fill this gap, in chapter 4 of this thesis, we develop a new paradigm to study a collective group of N agents moving and interacting in a three-dimensional space. Our paradigm helps to identify the spatio-temporal states of the motion of the group and their associated transition probabilities. This framework enables the estimation of the free energy landscape corresponding to the identified 75 states. Based on the energy landscape, we quantify missing information, emergence, self- organization and complexity for a collective motion. We show that the collective motion of the group of agents evolves to reach the most probable state with relatively lowest energy level and lowest missing information compared to other possible states. Our analysis demonstrates that the natural group of animals exhibit a higher degree of emergence, self-organization and complexity over time. Consequently, this algorithm can be integrated into new frameworks to engineer collective motions to achieve certain degrees of emergence, self-organization and complexity. • Recent Ebola outbreak in West Africa raised the need for developing a new treatment technique to control this deadly hemorrhagic virus. First step toward finding a cure for Ebola infection is to understand how the immune system interacts with the virus at the cellular level and determine robust control theoretic approaches for therapeutic purposes. In chapter 5 of this thesis, we propose a mathematical model for capturing the interactions among Ebola virions and immune cells and describe their dynamics through a set of coupled mean field equations. Our modeling framework provides a comprehensive perspective of the battle of Ebola virus against the host immune system. Based on our modeling formalism, we are able to study the effect of chemotherapy, gene therapy, and bacteria therapy on the network dynamics of virus and host immune system cells. This helps to design with the most efficient drug and minimize the side effects of trial-and-error approaches behind current drug design and therapeutic strategies. Finally, we demonstrate that bacteria therapy is a more reliable therapeutic strategy when compared to chemotherapy and gene therapy for Ebola treatment. 6.2 Future Research Direction Building on our work, we plan to propose a new mathematical model which captures the dynamics of a collective group of bacteria-propelled micro-robots. The bacteria-propelled motion enables these micro- robots to swim towards the inaccessible regions of human body for diagnosis and drug delivery purposes. The building blocks of such model are the mathematical characteristics of the bacterium motion, which we presented in chapter 2 of this thesis. More precisely, the model should follow anomalous diffusion and long range memory for the bacteria based on the nonergodicity and multi-fractality analysis. Moreover, based on the algorithm we developed in chapter 4, we plan to engineer the collective group of micro-robots to self-organize into required spatio-temporal arrangement based on the optimized energy landscape. Based on this, following our algorithmic strategy, one can set the parameters that governs the dynamic of one micro-robot and its corresponding interactions with other agents to achieve a certain degree of emergence, self-organization and complexity in an engineered micro-robot swarm. 76 7 Bibliography [1] W. C. Ruder, “Synthetic Biology Moving into the Clinic,” Science (80-. )., vol. 333, pp. 1248–1252, 2011. [2] O. Felfoul et al., “Magneto-aerotactic bacteria deliver drug-containing nanoliposomes to tumour hypoxic regions.,” Nat. Nanotechnol., no. August, 2016. [3] F. Matthäus, M. S. Mommer, T. Curk, and J. Dobnikar, “On the origin and characteristics of noise- induced Lévy walks of E. coli,” PLoS One, vol. 6, no. 4, p. e18623, 2011. [4] D. W. Stroock, “Some stochastic processes which arise from a model of the motion of a bacterium,” Zeitschrift f??r Wahrscheinlichkeitstheorie und Verwandte Gebiete, vol. 28, no. 4, pp. 305–315, 1974. [5] D. Selmeczi, S. Mosler, P. H. Hagedorn, N. B. Larsen, and H. Flyvbjerg, “Cell motility as persistent random motion: theories from experiments.,” Biophys. J., vol. 89, no. 2, pp. 912–931, 2005. [6] E. A. Codling, M. J. Plank, S. Benhamou, and J. R. S. Interface, “Random walk models in biology.,” J. R. Soc. Interface, vol. 5, no. 25, pp. 813–834, 2008. [7] T. Curk, D. Marenduzzo, and J. Dobnikar, “Chemotactic Sensing towards Ambient and Secreted Attractant Drives Collective Behaviour of E. coli,” PLoS One, vol. 8, no. 10, 2013. [8] A. Bricard, J.-B. Caussin, N. Desreumaux, O. Dauchot, and D. Bartolo, “Emergence of macroscopic directed motion in populations of motile colloids.,” Nature, vol. 503, no. 7474, pp. 95–98, 2013. [9] H. G. Othmer and C. Xue, “The Mathematical Analysis of Biological Aggregation and Dispersal: Progress, Problems and Perspectives,” Springer Berlin Heidelberg, 2013, pp. 79–127. [10] R. Erban and J. Haskovec, “From individual to collective behaviour of coupled velocity jump processes: A locust example,” Kinet. Relat. Model., vol. 5, no. 4, pp. 817–842, 2012. [11] C. Xue and H. Othmer, “Multiscale models of taxis-driven patterning in bacterial populations,” SIAM J. Appl. Math., vol. 70, no. 1, pp. 133–167, 2009. [12] T. Hillen and K. J. Painter, “A user’s guide to PDE models for chemotaxis,” Journal of Mathematical Biology, vol. 58, no. 1–2. pp. 183–217, 2009. [13] C. Escudero, “The fractional keller–segel model,” Nonlinearity, vol. 19, pp. 2909–2918, 2006. [14] A. Polezhaev, R. Pashkov, and A. Lobanov, “Spatial patterns formed by chemotactic bacteria Escherichia coli,” Int. J., 2003. 77 [15] M. E. Cates, “Diffusive transport without detailed balance in motile bacteria: does microbiology need statistical physics?,” Rep. Prog. Phys., vol. 75, no. 4, p. 42601, 2012. [16] M. Schnitzer, “Theory of continuum random walks and applications to chemotaxis,” Physical Review E, vol. 48, no. 4. pp. 2553–2568, 1993. [17] O. A. Croze, G. P. Ferguson, M. E. Cates, and W. C. K. Poon, “Migration of chemotactic bacteria in soft agar: Role of gel concentration,” Biophys. J., vol. 101, no. 3, pp. 525–534, 2011. [18] P. De Gennes, “Chemotaxis: the role of internal delays,” Eur. Biophys. J., vol. 33, pp. 691–693, 2004. [19] G. Ariel, A. Rabani, S. Benisty, J. D. Partridge, R. M. Harshey, and A. Be’er, “Swarming bacteria migrate by Lévy Walk.,” Nat. Commun., vol. 6, p. 8396, 2015. [20] J. Zhuang et al., “Analytical modeling and experimental characterization of chemotaxis in Serratia marcescens,” Phys. Rev. E, vol. 89, no. 5, p. 52704, 2014. [21] M. R. Edwards, R. W. Carlsen, J. Zhuang, and M. Sitti, “Swimming characterization of Serratia marcescens for bio-hybrid micro-robotics,” J. Micro-Bio Robot., vol. 9, no. 3, pp. 47–60, 2014. [22] U. Frisch, Turbulence: the legacy of AN Kolmogorov. Cambridge University Press, 1995. [23] P. Bogdan, M. Kas, R. Marculescu, and O. Mutlu, “QuaLe: A quantum-leap inspired model for non- stationary analysis of NoC traffic in chip multi-processors,” in NOCS 2010 - The 4th ACM/IEEE International Symposium on Networks-on-Chip, 2010, pp. 241–248. [24] N. Goldenfeld and C. Woese, “Life is Physics: Evolution as a Collective Phenomenon Far From Equilibrium,” Annu. Rev. Condens. Matter Phys., vol. 2, no. 1, pp. 375–399, 2011. [25] M. M. Dankulov, R. Melnik, and B. Tadić, “The dynamics of meaningful social interactions and the emergence of collective knowledge.,” Sci. Rep., vol. 5, no. July, p. 12197, 2015. [26] K. Tunstrøm, Y. Katz, C. Ioannou, and C. Huepe, “Collective states, multistability and transitional behavior in schooling fish,” PLoS Comput, vol. 9, no. 2, p. e1002915, 2013. [27] J. N. Taylor, C.-B. Li, D. R. Cooper, C. F. Landes, and T. Komatsuzaki, “Error-based extraction of states and energy landscapes from experimental single-molecule time-series.,” Sci. Rep., vol. 5, p. 9174, Mar. 2015. [28] I. Aoki, “A Simulation Study on the Schooling Mechanism in Fish,” Bull. Japanese Soc. Sci. Fish., vol. 48(8), no. 8, pp. 1081–1088, 1982. [29] C. W. Reynolds, “Flocks, herds and schools: A distributed behavioral model,” ACM SIGGRAPH Comput. Graph., vol. 21, no. 4, pp. 25–34, Aug. 1987. [30] Z. Csahk and T. Vicsek, “Lattice-gas model for collective biological motion,” Phys. Rev. E, vol. 52, no. 5, pp. 5297–5303, Nov. 1995. [31] S. Gueron, S. A. Levin, and D. I. Rubenstein, “The dynamics of herds: From individuals to 78 aggregations,” J. Theor. Biol., vol. 182, no. 1, pp. 85–98, 1996. [32] A. Czirók, M. Vicsek, and T. Vicsek, “Collective motion of organisms in three dimensions,” Phys. A Stat. Mech. its Appl., vol. 264, no. 1, pp. 299–304, 1999. [33] I. D. Couzin, J. Krause, N. R. Franks, and S. A. Levin, “Effective leadership and decision-making in animal groups on the move,” Nature, vol. 433, no. 7025, pp. 513–516, Feb. 2005. [34] I. Couzin, “Collective minds,” Nature, vol. 445, no. 7129, pp. 715–715, Feb. 2007. [35] T. Vicsek and A. Zafeiris, “Collective motion,” Physics Reports, vol. 517, no. 3–4. pp. 71–140, 2012. [36] R. Bale, M. Hao, A. P. S. Bhalla, and N. A. Patankar, “Energy efficiency and allometry of movement of swimming and flying animals.,” Proc. Natl. Acad. Sci. U. S. A., vol. 111, no. 21, pp. 7517–21, May 2014. [37] I. D. COUZIN, J. KRAUSE, R. JAMES, G. D. RUXTON, and N. R. FRANKS, “Collective Memory and Spatial Sorting in Animal Groups,” J. Theor. Biol., vol. 218, no. 1, pp. 1–11, 2002. [38] S. Hauert, J. C. Zufferey, and D. Floreano, “Evolved swarming without positioning information: An application in aerial communication relay,” Auton. Robots, vol. 26, no. 1, pp. 21–32, Jan. 2009. [39] S. Hauert, J. C. Zufferey, and D. Floreano, “Reverse-engineering of artificially evolved controllers for swarms of robots,” in 2009 IEEE Congress on Evolutionary Computation, CEC 2009, 2009, pp. 55–61. [40] B. T. Fine and D. A. Shell, “Eliciting collective behaviors through automatically generated environments,” in IEEE International Conference on Intelligent Robots and Systems, 2013, pp. 3303–3308. [41] D. A. Shell and M. J. Matarić, “On foraging strategies for large-scale multi-robot systems,” in IEEE International Conference on Intelligent Robots and Systems, 2006, pp. 2717–2723. [42] B. T. Fine and D. A. Shell, “Examining the Information Requirements for Flocking Motion,” Springer Berlin Heidelberg, 2012, pp. 442–452. [43] C. Nam and D. A. Shell, “When to do your own thing: Analysis of cost uncertainties in multi-robot task allocation at run-time,” 2015 IEEE Int. Conf. Robot. Autom., pp. 1249–1254, 2015. [44] N. AY and D. POLANI, “Information Flows in Causal Networks,” Adv. Complex Syst., vol. 11, no. 1, pp. 17–41, Feb. 2008. [45] J.-H. Kim and D. A. Shell, “A new model for self-organized robotic clustering: understanding boundary induced densities and cluster compactness,” Int. Conf. Robot. Autom., pp. 5858–5863, May 2015. [46] A. Martinoli, K. Easton, and W. Agassounon, “Modeling Swarm Robotic Systems: a Case Study in Collaborative Distributed Manipulation,” Int. J. Rob. Res., vol. 23, no. 4, pp. 415–436, 2004. 79 [47] K. Lerman, A. Galstyan, A. Martinoli, and A. Ijspeert, “A Macroscopic Analytical Model of Collaboration in Distributed Robotic Systems,” Artif. Life, vol. 7, no. 4, pp. 375–393, 2001. [48] L. Li, A. Martinoli, and Y. Abu-Mostafa, “Emergent specialization in swarm systems,” Intell. Data Eng. Autom. Learn., vol. 2412, pp. 261–266, 2002. [49] J. Halloy et al., “Social Integration of Robots into Groups of Cockroaches to Control Self-Organized Choices,” Science (80-. )., vol. 318, no. 5853, 2007. [50] G. Hummer and A. Szabo, “A Free energy reconstruction from nonequilibrium single-molecule pulling experiments,” Proc. Natl. Acad. Sci., vol. 98, no. 7, pp. 3658–3661, Mar. 2001. [51] P. Raiteri, A. Laio, F. L. Gervasio, C. Micheletti, and M. Parrinello, “Efficient reconstruction of complex free energy landscapes by multiple walkers metadynamics,” J. Phys. Chem. B, vol. 110, no. 8, pp. 3533–3539, 2006. [52] D. Prada-Gracia and J. Gómez-Gardeñes, “Exploring the free energy landscape: from dynamics to networks and back,” PLoS Comput, vol. 5, no. 6, p. e1000415, 2009. [53] P. Schuetz, R. Wuttke, B. Schuler, and A. Caflisch, “Free energy surfaces from single-distance information,” J. Phys. Chem. B, vol. 114, no. 46, pp. 15227–15235, Nov. 2010. [54] A. Baba and T. Komatsuzaki, “Extracting the underlying effective free energy landscape from single-molecule time series--local equilibrium states and their network.,” Phys. Chem. Chem. Phys., vol. 13, no. 4, pp. 1395–1406, 2011. [55] W. Jiang, Y. Luo, L. Maragliano, and B. Roux, “Calculation of free energy landscape in multi- dimensions with hamiltonian-exchange umbrella sampling on petascale supercomputer,” J. Chem. Theory Comput., vol. 8, no. 11, pp. 4672–4680, Nov. 2012. [56] D. J. Wales, “Decoding the energy landscape: extracting structure, dynamics and thermodynamics,” Philos. Trans. R. Soc. A Math. Phys. Eng. Sci., vol. 370, no. 1969, pp. 2877–2899, 2012. [57] S. Kawai and T. Komatsuzaki, “Effect of timescale on energy landscape: Distinction between free- energy landscape and potential of mean force,” Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., vol. 87, no. 3, p. 30803, Mar. 2013. [58] D. J. Wales and P. Salamon, “Observation time scale, free-energy landscapes, and molecular symmetry,” Proc. Natl. Acad. Sci., vol. 111, no. 2, pp. 617–622, Jan. 2014. [59] W. Li, W. Wang, and S. Takada, “Energy landscape views for interplays among folding, binding, and allostery of calmodulin domains.,” Proc. Natl. Acad. Sci. U. S. A., vol. 111, no. 29, pp. 10550– 5, Jul. 2014. [60] D. Mehta, C. Hughes, M. Kastner, and D. J. Wales, “Potential energy landscape of the two- dimensional XY model: Higher-index stationary points,” J. Chem. Phys., vol. 140, no. 22, p. 224503, Jun. 2014. 80 [61] L. C. Smeeton, M. T. Oakley, and R. L. Johnston, “Visualizing energy landscapes with metric disconnectivity graphs,” J. Comput. Chem., vol. 35, no. 20, pp. 1481–1490, Jul. 2014. [62] T. Watanabe, N. Masuda, F. Megumi, R. Kanai, and G. Rees, “Energy landscape and dynamics of brain activity during human bistable perception,” Nat. Commun., vol. 5, no. 1, p. 4765, Aug. 2014. [63] C. E. Shannon and C. E., “A mathematical theory of communication,” ACM SIGMOBILE Mob. Comput. Commun. Rev., vol. 5, no. 1, pp. 3–55, Jan. 2001. [64] V. Sourjik and H. C. Berg, “Receptor sensitivity in bacterial chemotaxis,” Pnas, vol. 99, no. 1, pp. 123–127, 2002. [65] C. Kim, M. Jackson, R. Lux, and S. Khan, “Determinants of chemotactic signal amplification in Escherichia coli.,” J. Mol. Biol., vol. 307, no. 1, pp. 119–135, 2001. [66] Y. Tu, “Quantitative modeling of bacterial chemotaxis: signal amplification and accurate adaptation.,” Annu. Rev. Biophys., vol. 42, no. February, pp. 337–359, 2013. [67] U. Alon, “An Introduction to Systems Biology: Design Principles of Biological Circuits,” Chapman HallCRC mathematical and computational biology series. CRC press, p. 301, 2006. [68] G. H. Wadhams and J. P. Armitage, “Making sense of it all: bacterial chemotaxis,” Nat.Rev.Mol.Cell Biol., vol. 5, no. 12, pp. 1024–1037, Dec. 2004. [69] P. Cluzel, M. Surette, and S. Leibler, “An ultrasensitive bacterial motor revealed by monitoring signaling proteins in single cells,” Science (80-. )., vol. 287, no. 5458, pp. 1652–1655, 2000. [70] M. M. McEvoy, A. Bren, M. Eisenbach, and F. W. Dahlquist, “Identification of the binding interfaces on CheY for two of its targets the phosphatase CheZ and the flagellar switch protein FliM,” J. Mol. Biol., vol. 289, no. 5, pp. 1423–1433, 1999. [71] M. Welch, K. Oosawa, S. Aizawa, and M. Eisenbach, “Phosphorylation-dependent binding of a signal molecule to the flagellar switch of bacteria.,” Proc. Natl. Acad. Sci. U. S. A., vol. 90, no. 19, pp. 8787–8791, 1993. [72] A. S. Toker and R. M. Macnab, “Distinct regions of bacterial flagellar switch protein FliM interact with FliG, FliN and CheY,” J. Mol. Biol., vol. 273, no. 3, pp. 623–634, 1997. [73] J. F. Hess, K. Oosawa, N. Kaplan, and M. I. Simon, “Phosphorylation of three proteins in the signalling pathway of bacterial chemotaxis,” Cell, vol. 53, no. 1, pp. 79–87, 1988. [74] G. S. Anand, P. N. Goudreau, and A. M. Stock, “Activation of methylesterase CheB: Evidence of a dual role for the regulatory domain,” Biochemistry, vol. 37, no. 40, pp. 14038–14047, Oct. 1998. [75] M. Li and G. L. Hazelbauer, “Cellular stoichiometries of the components of the chemotaxis signaling complex,” J. Bacteriol, vol. 186, no. 12, pp. 3687–3694, 2004. [76] J. R. Maddock and L. Shapiro, “Polar location of the chemoreceptor complex in the Escherichia coli cell.,” Science, vol. 259, no. 5102, pp. 1717–1723, 1993. 81 [77] J. E. Segall, S. M. Block, and H. C. Berg, “Temporal comparisons in bacterial chemotaxis.,” Proc. Natl. Acad. Sci. U. S. A., vol. 83, no. 23, pp. 8987–8991, 1986. [78] V. Sourjik and H. C. Berg, “Functional interactions between receptors in bacterial chemotaxis,” Nature, vol. 428, no. 6981, pp. 1–4, 2004. [79] D. Bray, “Bacterial chemotaxis and the question of gain.,” Proc. Natl. Acad. Sci. U. S. A., vol. 99, no. 1, pp. 7–9, 2002. [80] D. Bray, M. D. Levin, and C. J. Morton-Firth, “Receptor clustering as a cellular mechanism to control sensitivity.,” Nature, vol. 393, no. 6680, pp. 85–88, 1998. [81] R. Jasuja, J. Keyoung, G. P. Reid, D. R. Trentham, and S. Khan, “Chemotactic responses of Escherichia coli to small jumps of photoreleased L-aspartate.,” Biophys. J., vol. 76, no. 3, pp. 1706– 1719, 1999. [82] E. Ising, “Beitrag zur theorie des ferromagnetismus,” Zeitschrift f{ü}r Phys. A Hadron. Nucl., vol. 31, no. 1, pp. 253–258, 1925. [83] T. Duke and D. Bray, “Heightened sensitivity of a lattice of membrane receptors.,” Proc. Natl. Acad. Sci. U. S. A., vol. 96, no. 18, pp. 10104–10108, 1999. [84] B. a Mello and Y. Tu, “Perfect and near-perfect adaptation in a model of bacterial chemotaxis.,” Biophys. J., vol. 84, no. 5, pp. 2943–2956, 2003. [85] B. a Mello, L. Shaw, and Y. Tu, “Effects of receptor interaction in bacterial chemotaxis.,” Biophys. J., vol. 87, no. 3, pp. 1578–1595, 2004. [86] B. a Mello and Y. Tu, “An allosteric model for heterogeneous receptor complexes: understanding bacterial chemotaxis responses to multiple stimuli.,” Proc. Natl. Acad. Sci. U. S. A., vol. 102, no. 48, pp. 17354–17359, 2005. [87] B. a Mello and Y. Tu, “Effects of adaptation in maintaining high sensitivity over a wide range of backgrounds for Escherichia coli chemotaxis.,” Biophys. J., vol. 92, no. 7, pp. 2329–2337, 2007. [88] M. L. Skoge, R. G. Endres, and N. S. Wingreen, “Receptor-receptor coupling in bacterial chemotaxis: evidence for strongly coupled clusters.,” Biophys. J., vol. 90, no. 12, pp. 4317–4326, 2006. [89] G. Lan, S. Schulmeister, V. Sourjik, and Y. Tu, “Adapt locally and act globally: strategy to maintain high chemoreceptor sensitivity in complex environments.,” Mol. Syst. Biol., vol. 7, no. 1, p. 475, 2011. [90] J. Monod, J. Wyman, and J. P. Changeux, “on the Nature of Allosteric Transitions: a Plausible Model.,” J. Mol. Biol., vol. 12, no. 1, pp. 88–118, 1965. [91] J.-P. Changeux, “Allostery and the Monod-Wyman-Changeux model after 50 years.,” Annu. Rev. Biophys., vol. 41, pp. 103–133, 2012. 82 [92] C. V. Rao, M. Frenklach, and A. P. Arkin, “An allosteric model for transmembrane signaling in bacterial chemotaxis,” J. Mol. Biol., vol. 343, no. 2, pp. 291–303, 2004. [93] H. C. Berg and D. a Brown, “Chemotaxis in Escherichia coli analysed by three-dimensional tracking.,” Nature, vol. 239, no. 5374, pp. 500–504, Oct. 1972. [94] M. Kollmann, L. Lovdok, K. Bartholome, J. Timmer, and V. Sourjik, “Design principles of a bacterial signalling network,” Nature, vol. 438, no. 7067, pp. 504–507, 2005. [95] N. Barkai and S. Leibler, “Robustness in simple biochemical networks.,” Nature, vol. 387, no. 6636, pp. 913–917, 1997. [96] B. E. Knox, P. N. Devreotes, A. Goldbeter, and L. A. Segel, “A molecular mechanism for sensory adaptation based on ligand-induced receptor modification.,” Proc. Natl. Acad. Sci. U. S. A., vol. 83, no. 8, pp. 2345–2349, 1986. [97] L. Michaelis and M. L. Menten, “Die Kinetik der Invertinwirkung,” Biochem Z, vol. 49, no. February, pp. 333–369, 1913. [98] U. Alon, M. G. Surette, N. Barkai, and S. Leibler, “Robustness in bacterial chemotaxis.,” Nature, vol. 397, no. 6715, pp. 168–171, Jan. 1999. [99] T. M. Yi, Y. Huang, M. I. Simon, and J. Doyle, “Robust perfect adaptation in bacterial chemotaxis through integral feedback control.,” Proc. Natl. Acad. Sci., vol. 97, no. 9, pp. 4649–4653, 2000. [100] C. J. Morton-Firth, T. S. Shimizu, and D. Bray, “A free-energy-based stochastic simulation of the Tar receptor complex.,” J. Mol. Biol., vol. 286, no. 4, pp. 1059–1074, 1999. [101] Y. Meir, V. Jakovljevic, O. Oleksiuk, V. Sourjik, and N. S. Wingreen, “Precision and kinetics of adaptation in bacterial chemotaxis,” Biophys. J., vol. 99, no. 9, pp. 2766–2774, 2010. [102] Y. V Kalinin, L. J. L., Y. Tu, and M. Wu, “Logarithmic sensing in Escherichia coli bacterial chemotaxis,” Biophys. J., vol. 96, no. 6, pp. 2439–2448, 2009. [103] M. W. Sneddon, J. R. Faeder, and T. Emonet, “Efficient modeling, simulation and coarse-graining of biological complexity with NFsim,” Nat. Methods, vol. 8, no. 2, pp. 177–183, 2011. [104] T. Ahmed, T. S. Shimizu, and R. Stocker, “Bacterial chemotaxis in linear and nonlinear steady microfluidic gradients,” Nano Lett., vol. 10, no. 9, pp. 3379–3385, Sep. 2010. [105] N. Vladimirov, L. Løvdok, D. Lebiedz, and V. Sourjik, “Dependence of bacterial chemotaxis on gradient shape and adaptation rate,” PLoS Comput. Biol., vol. 4, no. 12, p. e1000242, Dec. 2008. [106] G. Wei, P. Bogdan, and R. Marculescu, “Efficient modeling and simulation of bacteria-based nanonetworks with bnsim,” IEEE J. Sel., vol. 31, no. 12, pp. 868–878, 2013. [107] Y. Tu, T. S. Shimizu, and H. C. Berg, “Modeling the chemotactic response of Escherichia coli to time-varying stimuli.,” Proc Natl Acad Sci U S A, vol. 105, no. 39, pp. 14855–14860, 2008. [108] A. Vaknin and H. C. Berg, “Physical Responses of Bacterial Chemoreceptors,” J. Mol. Biol., vol. 83 366, no. 5, pp. 1416–1423, Mar. 2007. [109] E. O. Budrene, Berg, and H. C. Berg, “Dynamics of formation of symmetrical patterns by chemotactic bacteria.,” Nature, vol. 376, no. 6535, pp. 49–53, 1995. [110] S.-Y. Cheng, S. Heilman, M. Wasserman, S. Archer, M. L. Shuler, and M. Wu, “A hydrogel-based microfluidic device for the studies of directed cell migration.,” Lab Chip, vol. 7, no. 6, pp. 763–769, 2007. [111] E. F. Keller and L. A. Segel, “Model for chemotaxis,” J. Theor. Biol., vol. 30, no. 2, pp. 225–234, 1971. [112] H. Park, C. C. Guet, T. Emonet, and P. Cluzel, “Fine-Tuning of Chemotactic Response in E. coli Determined by High-Throughput Capillary Assay,” Curr. Microbiol., vol. 62, no. 3, pp. 764–769, Mar. 2011. [113] K. C. Chen, R. M. Ford, and P. T. Cummings, “Mathematical models for motile bacterial transport in cylindrical tubes.,” J. Theor. Biol., vol. 195, no. 4, pp. 481–504, 1998. [114] I. Gomez Portillo, D. Campos, and V. Méndez, “Intermittent random walks: transport regimes and implications on search strategies,” J. Stat. Mech. Theory Exp., vol. 2011, no. 2, p. P02033, Feb. 2011. [115] F. Thiel, L. Schimansky-Geier, and I. M. Sokolov, “Anomalous diffusion in run-and-tumble motion,” Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., vol. 86, no. 2, p. 21117, Aug. 2012. [116] T. Ahmed and R. Stocker, “Experimental verification of the behavioral foundation of bacterial transport parameters using microfluidics.,” Biophys. J., vol. 95, no. 9, pp. 4481–4493, 2008. [117] H. Koorehdavoudi et al., “Multi-fractal characterization of bacterial swimming dynamics: a case study on real and simulated Serratia marcescens,” Proc. R. Soc. London A Math. Phys. Eng. Sci., vol. 473, no. 2203, 2017. [118] O. Chepizhko, E. G. Altmann, and F. Peruani, “Optimal noise maximizes collective motion in heterogeneous media,” Phys. Rev. Lett., vol. 110, no. 23, 2013. [119] O. Chepizhko and F. Peruani, “Diffusion, subdiffusion, and trapping of active particles in heterogeneous media,” Phys. Rev. Lett., vol. 111, no. 16, 2013. [120] L. Jiang, Q. Ouyang, and Y. Tu, “Quantitative Modeling of <italic>Escherichia coli</italic> Chemotactic Motion in Environments Varying in Space and Time,” PLoS Comput Biol, vol. 6, no. 4, p. e1000735, 2010. [121] S. B. Van Albada and P. R. Ten Wolde, “Differential affinity and catalytic activity of CheZ in E. coli chemotaxis,” PLoS Comput. Biol., vol. 5, p. e1000378, May 2009. [122] E. Korobkova, T. Emonet1, J. M. G. Vilar2, T. S. Shimizu, and P. Cluzel, “From molecular noise to behavioural variability in a single bacterium,” Nature, vol. 428, no. 6982, pp. 574–578, 2004. 84 [123] Y. Tu and G. Grinstein, “How White Noise Generates Power-Law Switching in Bacterial Flagellar Motors,” Phys. Rev. Lett., vol. 94, p. 208101, 2005. [124] J. Taktikos, V. Zaburdaev, and H. Stark, “Collective dynamics of model microorganisms with chemotactic signaling,” Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., vol. 85, no. 5, p. 51901, 2012. [125] A. Sengupta, S. Van Teeffelen, and H. Lowen, “Dynamics of a microorganism moving by chemotaxis in its own secretion,” Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., vol. 80, no. 3, pp. 0311221–0311229, 2009. [126] E. Lushi, R. E. Goldstein, and M. J. Shelley, “Auto-chemotactic active suspensions: modeling, analysis and simulations,” arXiv Prepr. arXiv1310.7614, 2013. [127] P. D’Ettorre and D. Hughes, Sociobiology of Communication: an Iterdisciplinar Perspective. Oxford University Press, 2008. [128] T. Emonet, C. M. Macal, M. J. North, C. E. Wickersham, and P. Cluzel, “AgentCell: A digital single- cell assay for bacterial chemotaxis,” Bioinformatics, vol. 21, no. 11, pp. 2714–2721, 2005. [129] B. Behkam and M. Sitti, “Bacterial flagella-based propulsion and on/off motion control of microscale objects,” Appl. Phys. Lett., vol. 90, no. 2, 2007. [130] N. Darnton, L. Turner, K. Breuer, and H. C. Berg, “Moving fluid with bacterial carpets.,” Biophys. J., vol. 86, no. 3, pp. 1863–1870, 2004. [131] M. S. Sakar et al., “Modeling, control and experimental characterization of microbiorobots,” Int. J. Rob. Res., vol. 30, no. 6, pp. 647–658, 2011. [132] M. R. Edwards, R. Wright Carlsen, and M. Sitti, “Near and far-wall effects on the three-dimensional motion of bacteria-driven microbeads,” Appl. Phys. Lett., vol. 102, no. 14, p. 143701, 2013. [133] C. Willert and M. Gharib, “Three-dimensional particle imaging with a single camera,” Exp. Fluids, vol. 12, no. 6, pp. 353–358, 1992. [134] M. Wu, J. W. Roberts, and M. Buckley, “Three-dimensional fluorescent particle tracking at micron- scale using a single camera,” Exp. Fluids, vol. 38, no. 4, pp. 461–465, 2005. [135] M. Wu, J. W. Roberts, S. Kim, D. L. Koch, and M. P. DeLisa, “Collective bacterial dynamics revealed using a three-dimensional population-scale defocused particle tracking technique,” Appl. Environ. Microbiol., vol. 72, no. 7, pp. 4987–4994, 2006. [136] S. Burov, J.-H. Jeon, R. Metzler, and E. Barkai, “Single particle tracking in systems showing anomalous diffusion: the role of weak ergodicity breaking.,” Phys. Chem. Chem. Phys., vol. 13, pp. 1800–1812, 2011. [137] C. Tsallis, S. V. F. Levy, A. M. C. Souza, and R. Maynard, “Statistical mechanics foundations of the ubiquity of {L}évy distribution in nature,” Phys. Rev. Lett., vol. 75, no. 20, pp. 3589–3592, 85 1995. [138] S. Burov, R. Metzler, and E. Barkai, “Aging and nonergodicity beyond the Khinchin theorem.,” Proc. Natl. Acad. Sci. U. S. A., vol. 107, no. 30, pp. 13228–13233, 2010. [139] A. Celani and M. Vergassola, “Bacterial strategies for chemotaxis response.,” Proc. Natl. Acad. Sci. U. S. A., vol. 107, no. 4, pp. 1391–6, 2010. [140] A. Celani, T. S. Shimizu, and M. Vergassola, “Molecular and Functional Aspects of Bacterial Chemotaxis,” J. Stat. Phys., vol. 144, no. 2, pp. 219–240, 2011. [141] K. J. Painter and T. Hillen, “Volume-filling and quorum-sensing in models for chemosensitive movement,” Can. Appl. Math. Q., vol. 10, no. 4, pp. 501–543, 2002. [142] R. Tyson, S. R. Lubkin, and J. D. Murray, “A minimal mechanism for bacterial pattern formation.,” Proc. Biol. Sci., vol. 266, no. 1416, pp. 299–304, 1999. [143] B. R. Parry, I. V. Surovtsev, M. T. Cabeen, C. S. O’Hern, E. R. Dufresne, and C. Jacobs-Wagner, “The bacterial cytoplasm has glass-like properties and is fluidized by metabolic activity,” Cell, vol. 156, no. 1–2, pp. 183–194, 2014. [144] R. Metzler, J.-H. Jeon, A. G. Cherstvy, and E. Barkai, “Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking.,” Phys. Chem. Chem. Phys., vol. 16, no. 44, pp. 24128–24164, 2014. [145] G. Bel and I. Nemenman, “Ergodic and non-ergodic anomalous diffusion in coupled stochastic processes,” New J. Phys., vol. 11, no. 8, p. 83009, Aug. 2009. [146] R. Klages, G. Radons, and I. Sokolov, Anomalous transport: foundations and applications. John Wiley & Sons, 2008. [147] F. Matthäus, M. Jagodić, and J. Dobnikar, “E. coli superdiffusion and chemotaxis - search strategy, precision and motility,” Biophys. J., vol. 97, no. 4, pp. 946–957, 2009. [148] I. M. Sokolov, “Models of anomalous diffusion in crowded environments,” Soft Matter, vol. 8, no. 35, pp. 9043–9052, 2012. [149] M. Javanainen et al., “Anomalous and normal diffusion of proteins and lipids in crowded lipid membranes,” Faraday Discuss., vol. 161, pp. 397–417, 2013. [150] E. Barkai, Y. Garini, and R. Metzler, “Strange kinetics of single molecules in living cells,” Phys. Today, vol. 65, no. 8, pp. 29–35, 2012. [151] A. A. Dubkov, B. Spagnolo, and V. V. Uchaikin, “LEVY FLIGHT SUPERDIFFUSION: AN INTRODUCTION,” Int. J. Bifurc. Chaos, vol. 18, no. 9, pp. 2649–2672, 2008. [152] G. Margolin and E. Barkai, “Nonergodicity of a Time Series Obeying L ´ evy Statistics,” J. Stat. Phys., vol. 122, no. 1, pp. 137–167, 2006. [153] R. Metzler and J. Klafter, “The random walk’s guide to anomalous diffusion: a fractional dynamics 86 approach,” Phys. Rep., vol. 339, no. 1, pp. 1–77, 2000. [154] K. V. Mardia, “Measures of multivariate skewness and kurtosis with applications,” Biometrika, vol. 57, no. 3, pp. 519–530, 1970. [155] K. V Mardia, “Applications of Some Measures of Multivariate Skewness and Kurtosis in Testing Normality and Robustness Studies,” Sankhyā Indian J. Stat. Ser. B, vol. 36, no. 2, pp. 115–128, 1974. [156] K. V Mardia and P. J. Zemroch, “Algorithm AS 84: Measures of multivariate skewness and kurtosis,” J. R. Stat. Soc. Ser. C (Applied Stat., vol. 24, no. 2, pp. 262–265, 1975. [157] N. Henze and B. Zirkler, “A class of invariant consistent tests for multivariate normality,” Commun. Stat. - Theory Methods, vol. 19, no. 10, pp. 3595–3617, 1990. [158] N. Henze and T. Wagner, “A New Approach to the BHEP Tests for Multivariate Normality,” J. Multivar. Anal., vol. 62, no. 1, pp. 1–23, 1997. [159] R. Johnson and D. Wichern, Applied multivariate statistical analysis. Englewood Cliffs, NJ: prentice Hall, 1992. [160] C. J. Mecklin and D. J. Mundfrom, “On Using Asymptotic Critical Values in Testing for Multivariate Normality,” InterStat, vol. 1, pp. 1–12, 2003. [161] C. J. Mecklin and D. J. Mundfrom, “A Monte Carlo comparison of the Type I and Type II error rates of tests of multivariate normality,” J. Stat. Comput. Simul., vol. 75, no. 2, pp. 93–107, 2005. [162] J. P. Royston, “An Extension of Shapiro and Wilk’s W Test for Normality to Large Samples,” J. R. Stat. Soc. Ser. C (Applied Stat., vol. 31, no. 2, pp. 115–124, 1982. [163] J. P. Royston, “Some techniques for assessing multivarate normality based on the Shapiro - Wilk W,” J. R. Stat. Soc., vol. 32, no. 2, pp. 121–133, 1983. [164] P. Royston, “Approximating the Shapiro-Wilk W-Test for non-normality,” Stat. Comput., vol. 2, no. 3, pp. 117–119, 1992. [165] P. Royston, “A Remark on Algorithm AS 181: The W-test for Normality,” J. R. Stat. Soc., vol. 44, no. 3, pp. 287–321, 1995. [166] M. B. Wilk, “An analysis of variance test for normality ( complete samples ) t,” Biometrika, vol. 52, no. 3, pp. 591–611, 2016. [167] R. B. D’agostino, “Transformation to normality of the null distribution of g1,” Biometrika, vol. 57, no. 3, pp. 679–681, 1970. [168] L. R. Shenton and K. O. Bowman, “A Bivariate Model for the Distribution of \sqrt{b}_1 and b_2,” J. Am. Stat. Assoc., vol. 72, no. 357, pp. 206–211, 1977. [169] J. A. Doornik and H. Hansen, “An Omnibus Test for Univariate and Multivariate Normality,” Oxford Bull. Econ., pp. 1–16, 1994. 87 [170] T. Svantesson and J. Wallace, “Tests for assessing multivariate normality and the covariance structure of MIMO data,” 2003 IEEE Int. Conf. Acoust. Speech, Signal Process. 2003. Proceedings. (ICASSP ’03)., vol. 4, pp. 656–659, 2003. [171] E. A. F. Ihlen, “Introduction to multifractal detrended fluctuation analysis in Matlab,” Front. Physiol., vol. 3 JUN, 2012. [172] J. W. Kantelhardt, S. A. Zschiegner, E. Koscielny-Bunde, S. Havlin, A. Bunde, and H. E. Stanley, “Multifractal detrended fluctuation analysis of nonstationary time series,” Physica A, vol. 316, no. 1, pp. 87–114, 2002. [173] B. Mandelbrot, Les Objects Fractals: Forme, Hasard et Dimention. Flammarion, Paris, 1984. [174] S. Timashev, “Self-similarity in Nature,” Stoch. CHAOTIC Dyn., vol. 502, pp. 562–566, 2000. [175] K. J. Falconer, Fractal geometry : mathematical foundations and applications. John Wiley & Sons, 2003. [176] A. Eke, P. Herman, B. G. Sanganahalli, F. Hyder, P. Mukli, and Z. Nagy, “Pitfalls in fractal time series analysis : fMRI BOLD as an exemplary case Pitfalls in fractal time series analysis : fMRI BOLD as an exemplary case,” Front., vol. 3, no. 417, p. 73, 2012. [177] A. Sokolov and I. S. Aranson, “Reduction of viscosity in suspension of swimming bacteria,” Phys. Rev. Lett., vol. 103, no. 14, 2009. [178] E. Ben-Jacob, O. Schochet, A. Tenenbaum, I. Cohen, A. Czirók, and T. Vicsek, “Generic modelling of cooperative growth patterns in bacterial colonies.,” Nature, vol. 368, no. 6466, pp. 46–49, Mar. 1994. [179] E. Ben-Jacob and H. Levine, “Self-engineering capabilities of bacteria.,” J. R. Soc. Interface, vol. 3, no. 6, pp. 197–214, 2006. [180] P. Bogdan, B. M. Deasy, B. Gharaibeh, T. Roehrs, and R. Marculescu, “Heterogeneous structure of stem cells dynamics: statistical models and quantitative predictions.,” Sci. Rep., vol. 4, p. 4826, 2014. [181] U. Mitra, N. Michelusi, S. Pirbadian, H. Koorehdavoudi, M. Y. El-Naggar, and P. Bogdan, “Queueing theory as a modeling tool for bacterial interaction: Implications for microbial fuel cells,” in 2015 International Conference on Computing, Networking and Communications, ICNC 2015, 2015, pp. 658–662. [182] P. Bogdan, “Mathematical Modeling and Control of Multifractal Workloads for Data-Center-on-a- Chip Optimization,” Proc. 9th Int. Symp. Networks-on-Chip, p. 21:1--21:8, 2015. [183] H. Koorehdavoudi and P. Bogdan, “A Statistical Physics Characterization of the Complex Systems Dynamics: Quantifying Complexity from Spatio-Temporal Interactions,” Sci. Rep., vol. 6, no. January, p. 27602, 2016. 88 [184] D. Chu, R. Strand, and R. Fjelland, “Theories of complexity: Common denominators of complex systems,” Complexity, vol. 8, no. 3. Wiley Subscription Services, Inc., A Wiley Company, pp. 19– 30, Jan-2003. [185] J. Jost, “External and internal complexity of complex adaptive systems,” Theory in Biosciences, vol. 123, no. 1. Springer-Verlag, pp. 69–88, Jun-2004. [186] C. B. Li, H. Yang, and T. Komatsuzaki, “New quantification of local transition heterogeneity of multiscale complex networks constructed from single-molecule time series,” J. Phys. Chem. B, vol. 113, no. 44, pp. 14732–14741, Nov. 2009. [187] C. Adami, “What is complexity?,” BioEssays, vol. 24, no. 12. Wiley Subscription Services, Inc., A Wiley Company, pp. 1085–1094, Dec-2002. [188] S. Chen et al., “Computing Free Energy Landscapes : Application to Ni-based Electrocatalysts with Pendant Amines for H2 Production and Oxidation Computing Free Energy Landscapes : Application to Ni-based Electrocatalysts with Pendant Amines for H 2 Production and Oxidatio,” ACS, vol. 4, no. 1, pp. 229–242, 2013. [189] C. A. Floudas and P. M. Pardalos, Optimization in computational chemistry and molecular biology : local and global approaches, no. 40. Springer Science & Business Media, 2000. [190] V. . Arnold, “Mathematical methods of classical mechanics,” Springer-Verlag, vol. 60. Springer New York, New York, NY, pp. 1–536, 1989. [191] H. Eyring, “The Energy of Activation for Bimolecular Reactions Involving Hydrogen and the Halogens, According to the Quantum Mechanics,” J. Am. Chem. Soc., vol. 53, no. 7, pp. 2537–2549, 1931. [192] K. a Dill, “Theory for the folding and stability of globular proteins.,” Biochemistry, vol. 24, no. APRIL 1985, pp. 1501–1509, 1985. [193] A. Fernández and E. I. Shakhnovich, “Activation-energy landscape for metastable RNA folding,” Phys. Rev. A, vol. 42, no. 6, pp. 3657–3659, 1990. [194] M. Mann, M. Kucharík, C. Flamm, and M. T. Wolfinger, “Memory efficient RNA energy landscape exploration.,” Bioinformatics, vol. 30, no. 18, pp. 2584–2591, 2014. [195] H. S. Seung, “How the brain keeps the eyes still.,” Proc. Natl. Acad. Sci. U. S. A., vol. 93, no. 23, pp. 13339–13344, 1996. [196] A. Baba and T. Komatsuzaki, “Construction of effective free energy landscape from single-molecule time series.,” Proc. Natl. Acad. Sci. U. S. A., vol. 104, no. 49, pp. 19297–19302, Dec. 2007. [197] L. V. Kantorovich and V. I. Kyrlov, Approximate method of higher analysis. 1958. [198] C. Gershenson and N. Fernández, “Complexity and information: Measuring emergence, self- organization, and homeostasis at multiple scales,” Complexity, vol. 18, no. 2, pp. 29–44, 2012. 89 [199] M. Gell-Mann and S. Lloyd, “Information measures, effective complexity, and total information,” Complexity, vol. 2, no. 1, pp. 44–52, 1996. [200] J. M. McGloin, C. J. Sullivan, and L. W. Kennedy, When Crime Appears: The Role of Emergence. Routledge, 2011. [201] A. Berdahl, C. J. Torney, C. C. Ioannou, J. J. Faria, and I. D. Couzin, “Emergent Sensing of Complex Environments by Mobile Animal Groups,” Science (80-. )., vol. 339, no. 6119, pp. 574–576, 2013. [202] H. Feldmann et al., “Effective post-exposure treatment of ebola infection,” PLoS Pathog., vol. 3, no. 1, pp. 0054–0061, 2007. [203] Y. Bar-Yam, “Multiscale Complexity/Entropy,” Adv. Complex Syst., vol. 7, no. 1, pp. 47–63, Mar. 2004. [204] M. Nagy et al., “Hierarchical group dynamics in pigeon flocks.,” Nature, vol. 464, no. 7290, pp. 890–893, Apr. 2010. [205] T. A. O’Shea-Wheller, A. B. Sendova-Franks, and N. R. Franks, “Differentiated anti-predation responses in a superorganism,” PLoS One, vol. 10, no. 11, p. e0141012, Nov. 2015. [206] B. Hölldobler and E. O. Wilson, “The Superorganism: The Beauty, Elegance, and Strangeness of Insect Societies,” Nature, vol. 456, no. January, p. 544, 2009. [207] F. Poiesi and A. Cavallaro, “Tracking multiple high-density homogeneous targets,” IEEE Trans. Circuits Syst. Video Technol., vol. 25, no. 4, pp. 623–637, 2015. [208] S. B. Rosenthal, C. R. Twomey, A. T. Hartnett, H. S. Wu, and I. D. Couzin, “Revealing the hidden networks of interaction in mobile animal groups allows prediction of complex behavioral contagion,” Proc. Natl. Acad. Sci., vol. 112, no. 15, pp. 4690–4695, Apr. 2015. [209] G. Beni and G. Beni, “From Swarm Intelligence to Swarm Robotics,” Robotics, vol. 3342, pp. 1–9, 2005. [210] J. Krause and G. D. Ruxton, “Living in Groups,” Living in groups, vol. I, no. 4. p. 210, 2002. [211] N. Tishby, F. Pereira, and W. Bialek, “The information bottleneck method,” arXiv Prepr. physics/0004057, pp. 368–377, 1999. [212] G. Ausiello, P. Crescenzi, G. Gambosi, and V. Kann, Complexity and approximation: Combinatorial optimization problems and their approximability properties. Springer Science & Business Media, 2012. [213] L. Fortnow, “The status of the P versus NP problem,” Commun. ACM, vol. 52, pp. 78–86, 2009. [214] K. Lerman, A. Martinoli, and A. Galstyan, “A Review of Probabilistic Macroscopic Models for Swarm Robotic Systems,” Swarm Robot., vol. 3342, no. i, pp. 143–152, 2004. [215] J. S. Weitz and J. Dushoff, “Modeling post-death transmission of Ebola: challenges for inference and opportunities for control.,” Sci. Rep., vol. 5, p. 8751, 2015. 90 [216] C. Peters and J. Peters, “An introduction to Ebola: the virus and the disease,” J. Infect. Dis., vol. 179, 1999. [217] D. Kirschner, “Using mathematics to understand HIV immune dynamics,” AMS Not., vol. 43, no. 2, pp. 191–202, 1996. [218] A. Yonezawa, M. Cavrois, and W. C. Greene, “Studies of ebola virus glycoprotein-mediated entry and fusion by using pseudotyped human immunodeficiency virus type 1 virions: involvement of cytoskeletal proteins and enhancement by tumor necrosis factor alpha.,” J. Virol., vol. 79, no. 2, pp. 918–926, 2005. [219] H. C. Tuckwell, Tuckwell, and H. C., “Viral Population Growth Models,” in Encyclopedia of Biostatistics, Chichester, UK: John Wiley & Sons, Ltd, 2005. [220] S. Banton, Z. Roth, and M. Pavlovic, “Mathematical Modeling of Ebola Virus Dynamics as a Step towards Rational Vaccine Design.pdf,” IFMBE Proceedings. pp. 196–200, 2010. [221] D. G. Bausch, A. G. Sprecher, B. Jeffs, and P. Boumandouki, “Treatment of Marburg and Ebola hemorrhagic fevers: A strategy for testing new drugs and vaccines under outbreak conditions,” Antiviral Res., vol. 78, no. 1, pp. 150–161, Apr. 2008. [222] D. Kirschner, S. Lenhart, and S. Serbin, “Optimal control of the chemotherapy of HIV.,” J. Math. Biol., vol. 35, no. 7, pp. 775–792, 1997. [223] P. Ankomah and B. R. Levin, “Two-drug antimicrobial chemotherapy: A mathematical model and experiments with Mycobacterium marinum,” PLoS Pathog., vol. 8, no. 1, p. e1002487, Jan. 2012. [224] E. De Clercq, “New developments in anti-HIV chemotherapy,” Biochim. Biophys. Acta (BBA)- Molecular Basis, vol. 1587, no. 2, pp. 258–275, 2002. [225] S. Aviran, P. S. Shah, D. V. Schaffer, and A. P. Arkin, “Computational models of HIV-1 resistance to gene therapy elucidate therapy design principles,” PLoS Comput. Biol., vol. 6, no. 8, p. e1000883, Aug. 2010. [226] D. L. DiGiusto et al., “RNA-based gene therapy for HIV with lentiviral vector-modified CD34(+) cells in patients undergoing transplantation for AIDS-related lymphoma.,” Sci. Transl. Med., vol. 2, no. 36, p. 36ra43, 2010. [227] S. T. Rutherford and B. L. Bassler, “Bacterial quorum sensing: its role in virulence and possibilities for its control.,” Cold Spring Harbor perspectives in medicine, vol. 2, no. 11. p. a012427, 2012. [228] O. Martinez et al., “Zaire Ebola virus entry into human dendritic cells is insensitive to cathepsin L inhibition,” Cell Microbiol, vol. 12, no. 2, pp. 148–157, 2010. [229] G. Ruthel et al., “Association of ebola virus matrix protein VP40 with microtubules.,” J. Virol., vol. 79, no. 8, pp. 4709–4719, 2005. [230] D. Dube et al., “Cell adhesion promotes Ebola virus envelope glycoprotein-mediated binding and 91 infection,” J Virol, vol. 82, no. 14, pp. 7238–7242, 2008. [231] J. Pettitt et al., “Therapeutic intervention of Ebola virus infection in rhesus macaques with the MB- 003 monoclonal antibody cocktail,” Sci Transl Med, vol. 5, no. 199, p. 199ra113, 2013. [232] T. Maruyama et al., “Ebola virus can be effectively neutralized by antibody produced in natural human infection.,” J. Virol., vol. 73, no. 7, pp. 6024–6030, 1999. [233] D. M. Tscherne, B. Manicassamy, and A. García-Sastre, “An enzymatic virus-like particle assay for sensitive detection of virus entry,” J. Virol. Methods, vol. 163, no. 2, pp. 336–343, 2010. [234] P. C. May et al., “The role of cathepsins in ocular physiology and pathology,” Exp Eye Res, vol. 84, no. 3, pp. 383–388, 2007. [235] T. S. Shimizu, N. Delalez, K. Pichler, and H. C. Berg, “Monitoring bacterial chemotaxis by using bioluminescence resonance energy transfer: absence of feedback from the flagellar motors.,” Proc. Natl. Acad. Sci. U. S. A., vol. 103, no. 7, pp. 2093–2097, 2006.

Abstract (if available)

Abstract Microbial communities and their interactions within the community and across communities and their environment is fundamental for evolution and survival of multicellular systems. For example, the heterogeneous types and interactions between gut bacteria are essential for digestion, extracting energy from raw material and even influencing the emergence and progression of diseases or the variations in psychological behavior (e.g., influencing the mood of a person). In addition, scientists showed that some bacteria can power electronic devices and build complex nanostructures. In medicine, one can exploit the intelligent bacteria motility and engineer swarms in order to detect different diseases and deliver drugs. To harness the intelligence and capabilities of microbial communities, we first have to understand their complex dynamics and be able to predict their motion in an interactive environment. This helps to improve the modeling of their motion at the cell-level and population-level. ❧ Traditional modeling of bacteria populations relies on simplifying assumptions such as the memoryless dynamics of their motion and ignores most of the physical and chemical interactions among them. To improve this state of affairs, we study the dynamics of bacterial communities in vitro and in silico, while accounting for realistic conditions like volume exclusion, chemical interactions, obstacles, and distribution of chemoattractant in the environment. We demonstrate that the bacterial swarm dynamics is not only nonergodic but also nonlinear and multi-fractal in nature. Finally, we demonstrate that bacteria dynamics behavior strongly depends on the density and chemoattractant concentration. These mathematical characteristics represent the canonical set of features that need to be considered for the accurate modeling of bacteria motion based on the nature of their swarm dynamics. ❧ To quantify the degree of complexity in microbial communities, we develop a new paradigm to study a collective group of N agents moving and interacting in a three-dimensional space. Our paradigm helps to identify the spatio-temporal states of the motion of the group and their associated transition probabilities. This framework enables the estimation of the free energy landscape corresponding to the identified states. Based on the energy landscape, we quantify missing information, emergence, self-organization and complexity for a collective motion. We show that the collective motion of the group of agents evolves to reach the most probable state with relatively lowest energy level and lowest missing information compared to other possible states. Our analysis demonstrates that the natural group of animals exhibit a higher degree of emergence, self-organization and complexity over time. Consequently, this algorithm can be integrated into new frameworks to engineer bacteria propelled micro-robots to achieve certain degrees of emergence, self-organization and complexity and perform drug-delivery task inside human body. ❧ In the end, we propose bacteria therapy to cure Ebola infection as an application of using bacteria for therapeutic purposes. We developed mathematical model of interacting Ebola virus with human immune cells during bacteria therapy. Overall, our new contributions can help to engineer a swarm of bacteria for drug delivery purposes.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Microbial interaction networks: from single cells to collective behavior

PDF

A meta-interaction model for designing cellular self-organizing systems

PDF

Transitions in the collective behavior of microswimmers confined to microfluidic chips

PDF

Theoretical and computational foundations for cyber‐physical systems design

PDF

Investigations of fuel and hydrodynamic effects in highly turbulent premixed jet flames

PDF

Understanding dynamics of cyber-physical systems: mathematical models, control algorithms and hardware incarnations

PDF

Transient modeling, dynamic analysis, and feedback control of the Inductrack Maglev system

PDF

Heterogeneous graphs versus multimodal content: modeling, mining, and analysis of social network data

PDF

Modeling and dynamic analysis of coupled structure-moving subsystem problem

PDF

Using nonlinear feedback control to model human landing mechanics

PDF

Modeling and analysis of parallel and spatially-evolving wall-bounded shear flows

PDF

New approaches in modeling and control of dynamical systems

PDF

Dynamic social structuring in cellular self-organizing systems

PDF

Advancing energy recovery from food waste using anaerobic biotechnologies: performance and microbial ecology

PDF

Understanding human-building interactions through perceptual decision-making processes

PDF

Parametric and non‐parametric modeling of autonomous physiologic systems: applications and multi‐scale modeling of sepsis

PDF

Mechanical and flow interactions facilitate cooperative transport and collective locomotion in animal groups

PDF

Feature and model based biomedical system characterization of cancer

PDF

Enabling human-building communication to promote pro-environmental behavior in office buildings

PDF

Computational tumor ecology: a model of tumor evolution, heterogeneity, and chemotherapeutic response

Asset Metadata

Creator Koorehdavoudi, Hana (author)

Core Title Mathematical characterizations of microbial communities: analysis and implications

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Mechanical Engineering

Publication Date 07/26/2017

Defense Date 12/01/2016

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag energy landscape,mathematical characterizations,microbial communities,multi-fractality,OAI-PMH Harvest,single bacterium motion,swarm motion

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Flashner, Henryk (committee chair), Bogdan, Paul (committee member), Spedding, Geoff (committee member)

Creator Email h.koorehdavoudi@gmail.com,koorehda@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-418642

Unique identifier UC11213612

Identifier etd-Koorehdavo-5646.pdf (filename),usctheses-c40-418642 (legacy record id)

Legacy Identifier etd-Koorehdavo-5646.pdf

Dmrecord 418642

Document Type Dissertation

Format application/pdf (imt)

Rights Koorehdavoudi, Hana

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA