An agent-based model of Trypanosoma brucei social motility to explore determinants of colony pattern formation

An agent-based model of Trypanosoma brucei social motility to explore determinants of colony pattern formation | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results An agent-based model of Trypanosoma brucei social motility to explore determinants of colony pattern formation View ORCID Profile Andreas Kuhn , Timothy Krüger , Markus Engstler , View ORCID Profile Sabine C. Fischer doi: https://doi.org/10.1101/2025.11.04.686461 Andreas Kuhn 1 Julius-Maximilians-Universität Würzburg, Biocenter, Chair for Computational and Theoretical Biology (CCTB) , Klara-Oppenheimer-Weg 32, 97074 Würzburg, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Andreas Kuhn Timothy Krüger 2 Julius-Maximilians-Universität Würzburg, Biocenter, Department of Cell and Developmental Biology , Am Hubland, 97074 Würzburg, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site Markus Engstler 2 Julius-Maximilians-Universität Würzburg, Biocenter, Department of Cell and Developmental Biology , Am Hubland, 97074 Würzburg, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site Sabine C. Fischer 1 Julius-Maximilians-Universität Würzburg, Biocenter, Chair for Computational and Theoretical Biology (CCTB) , Klara-Oppenheimer-Weg 32, 97074 Würzburg, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Sabine C. Fischer For correspondence: sabine.fischer{at}uni-wuerzburg.de Abstract Full Text Info/History Metrics Preview PDF Abstract In vitro colonies of the unicellular parasite Trypanosoma brucei expand radially and establish fingering instabilities, a collective behavior known as social motility. The underlying mechanisms are thought to involve single-cell motility, chemical communication among cells, and mechanical interactions with the liquid boundary, but their relative contributions remain unclear. We aimed to determine which of the mechanisms are necessary to quantitatively reproduce the morphological characteristics of social motility. We developed a two-dimensional agent-based model that simulates colonies of 10 5 − 10 6 cells at single-cell resolution—two to four orders of magnitude larger than previous models. Cells are represented as point particles executing directional random walks with auto-chemotactic alignment, combined with reflective and mutable boundaries and exponential colony growth. The model was quantitatively evaluated by applying our previously established morphology metrics. We show quantitative agreement of the simulation results and experimental data in terms of colony morphology. Parameter exploration revealed that finger formation arises within a narrow range of trypanosome motility parameters that balance stochasticity and alignment, while boundary conditions modulate the speed of colony expansion. The diffusion coefficient of the chemotactic signal is the key determinant of pattern formation. Realistic behavior occurs at 2× 10 −11 – 10 −10 m 2 /s which corresponds to molecules of 12.1–1690 kDa. These results demonstrate that complex colony morphologies can emerge from minimal cell-level rules, suggesting testable hypotheses for the molecular drivers of trypanosome social motility. Furthermore, our approach provides a framework for dissecting the interplay between motility, signaling, and mechanical confinement in other microbial systems exhibiting collective behavior. Author summary In vitro colonies of the unicellular parasite Trypanosoma brucei expand radially and form finger-like patterns, a behavior known as social motility. The ability of trypanosomes to exhibit social motility has been linked to their successful journey through their insect vector, the tsetse fly, which facilitates their survival and spreading between hosts. However, the mechanisms driving this behavior are not yet clear. Experimental data suggest that the movement and growth of individual cells, their interactions with each other, and their effects on colony boundaries may all play a role. These different factors are difficult to separate in experiments. We employed mathematical modeling to investigate the relative importance of these factors and developed a model that represents individual cells. The simulations reproduced finger formation patterns similar to those observed in experiments. The results show that complex colony shapes can emerge from simple cell behaviors. We also found that the speed at which a signaling chemical spreads between cells is crucial, and that the predicted values do not match the properties of the chemicals that have been proposed so far. Identifying and testing candidate signaling molecules in experiments could be the next step. Additionally, our approach may also aid in understanding collective behaviors in other microorganisms. Introduction Understanding how simple individual behaviour gives rise to complex collective phenomena is a fundamental challenge across biology. Social motility (SoMo) [ 1 – 3 ] in trypanosomes presents a striking example of this emergence: densely packed individual cells, swimming in an apparently random way [ 4 ] on agarose gel surfaces, form intricate colony-level patterns. Initially, these colonies exhibit a radial outgrowth and, in a second phase, produce characteristic fingering instabilities [ 5 ]. The fingers or protrusions grow in a distinct perpendicular alignment away from the colony’s centre. While the molecular and cellular mechanisms underlying individual trypanosome motility have been extensively studied [ 6 – 9 ], their collective behaviour is still under investigation. Current hypotheses propose that social motility results from chemotactic responses at the colony level [ 10 – 12 ], with cells reacting to environmental cues such as exosomes [ 10 ], neighbouring E. coli colonies [ 11 ], or self-generated pH gradients resulting from glucose metabolism [ 12 , 13 ]. Alternative hypotheses suggest that physical interactions, including hydrodynamic coupling, steric alignment between adjacent cells, or interactions with the liquid–solid boundary of the colony, contribute to the emergent patterns observed during collective movement [ 4 ]. The relative contributions of these chemical and physical mechanisms are unclear. Addressing this question requires an individual-based model. Despite visual similarities between T. brucei colony morphologies and those observed in bacterial swarming systems [ 14 – 20 ], it remains unclear to what extent insights from bacterial motility can be transferred. Such bacterial systems have been modelled using continuum approaches, in which population-level behaviour is described through density fields governed by reaction–diffusion equations [ 16 , 17 , 21 , 22 ] or by individual-based models [ 18 , 23 ]. Pseudomonas aeruginosa stands out as the system most closely resembling T. brucei in both colony morphology with similar finger instabilities [ 24 ] and as well in microswimmer dynamics with similar swimming modes [ 25 , 26 ] as trypoansoma brucei [ 27 , 28 ]. However, a mechanistic understanding of P. aeruginosa swarming remains incomplete and is complicated by system-specific features, such as the secretion of rhamnolipids, which change the physical properties of the colony boundary. This limits the direct transferability of concepts to T. brucei . These limitations indicate the requirement for a new individual-based modelling approach to identify the microscopic mechanisms that drive social motility in trypanosome colonies. The two-dimensional agent-based model simulates individual trypanosomes at populations of ∼ 10 5 to 10 6 cells, consistent with social motility assays [ 1 , 2 , 4 , 5 , 10 ]. This represents an increase of four orders of magnitude in agent numbers compared to existing models of trypanosome [ 27 – 29 ], which examine the swimming dynamics of individual cells, and an increase of two orders of magnitude compared to other microswimmer models [ 30 – 33 ] for collective behaviour. We employ several key assumptions to make such large-scale simulations computationally feasible over the social motility timescales of up to 3 days. First, we represent each trypanosome as a point particle on a square lattice, with spacing matched to the projected area of real trypanosomes. The complex swimming motion of individual trypanosomes is reduced to a dry active matter model [ 34 ] where agents perform a directional random walk with reflective boundary conditions at the colony boundaries, capturing the essential persistent motion of these microorganisms in confined spaces [ 11 ] while abstracting away their detailed swimming dynamics [ 35 , 36 ]. Second, we model the complex liquid boundary of the colony on top of the semi-liquid agarose matrix by dividing the colony space into walkable and non-walkable lattice sites. Collisions of the agents with the non-walkable lattice sites on the colony boundary make those sites progressively walkable and part of the colony. Such an approach has been successfully used before to model the growth of bacterial colonies [ 22 , 37 ]. Third, we model the proposed chemical signalling and diffusion that mediates cell-cell communication by a reaction-diffusion equation and negative autochemotaxis similar to previous models for bacteria chemotaxis [ 16 , 21 , 37 – 40 ]. For this, we employ a separate, coarser grid. Our dual-grid approach balances biological realism with computational efficiency, allowing us to capture diffusion-mediated interactions while maintaining tractable simulation times. Through analysis and comparison with experimental observations using our previously established metrics to quantify Trypanosoma brucei colonies [ 5 ], we demonstrate that these simplified rules enable the reproduction of the morphological characteristics of social motility. We further analyse how the individual components change the emerging behaviour and found that they exhibit distinct effects. The colony morphology is very sensitive to the motility parameters of the single cells. The parameters for boundary interactions mostly affect expansion speed. The diffusion coefficient of the chemotactic signal affects both expansion speed and colony morphology. Hence, it is a key parameter for obtaining finger-like patterns comparable to experiments. We could not identify the chemical component that most likely drives Social motility because the optimal diffusion parameter value does not match any of the previously suggested chemical components. Materials and methods Model Architecture To model trypanosome social motility, we combined an agent-based model for the individual cells with a continuous model for chemical signalling and diffusion. The state of the model is stored in three primary data structures ( Figure 1 ). First, a table of agents representing the trypanosomes contains for each entry a unique identifier (id), the position x ( t ), and the orientation θ ( t ) as the angle in the 2D plane relative to the horizontal axis ( Figure 1 (a) ). Second, a 2D integer array called the agent-based model grid (ABM grid) models the spatial environment. Third, a 2D floating-point array called the gradient grid stores the chemical concentration values s ( x ). The ABM grid is a square lattice for which the lattice constant Δ x was set to 3.75 µ m, corresponding to a unit cell area of 14.06 µ m 2 (see Section S.1 for more details). Each grid cell stores an integer value k , where: Download figure Open in new tab Fig 1. Schematic of the model architecture: (a) Conceptual representation of individual components: ABM grid (left), agent table/list (right), and gradient grid ( s f = 2) (below). Each agent (black dot) occupies one cell in the ABM grid; walkable cells are white, and non-walkable cells are green. The agent table stores the position x ( t ) and direction θ ( t ) (blue arrow) of all agents. With scaling factor s f = 2, four ABM grid cells map to one gradient grid cell, with chemical concentration indicated by blue colour intensity. (b) Computational implementation of the data structures as in (a): The ABM grid is a 6×6 2D integer array, the agent table is a vector with position (integer tuple) and orientation (floating-point between 0° and 360°), and the gradient grid is a smaller 3×3 2D floating-point array. (c) Conceptual representation of integrated model: agents move in their current direction through varying chemical concentrations on empty ABM grid cells. k > 0 indicates the unique ID of an agent occupying that cell (only one agent per grid cell is possible), k = 0 represents an empty, walkable grid cell, k < 0 denotes a non-walkable grid cell outside the colony, with the absolute value of k determining the grid strength. The gradient grid stores the local chemical concentration at each grid point and handles calculations for chemical diffusion, adsorption, and decay. This grid can be scaled down by an integer scaling factor s f (e.g., 1, 2, 3, 4, …) relative to the ABM grid. Therefore, each gradient grid cell corresponds to ABM grid cells. Initialization The model initialisation proceeds in several sequential steps: Creation of the ABM grid with desired size L and corresponding number of grid points. Assignment of the same negative value k 0 (grid strength) to all grid points. Definition of the initial colony geometry by setting all grid points to zero within a circle with radius r colony away from the center of the ABM grid. Random placement of N agents on walkable grid points ( k = 0) with random orientations drawn from a uniform distribution. Creation of the gradient grid, scaled to the ABM grid with factor s f and initialized with a small, randomly distributed chemical concentration at each grid point. Simulations The simulations are performed in an iterative manner. Updates of the gradient grid, the agents, and the ABM grid are combined ( Figure 2 ). In the following, we describe the individual parts in detail. Download figure Open in new tab Fig 2. Flowchart of the internal model logic. The blue ovals represent functions that are called in consecutive order, and the yellow boxes represent the objects that they modify. The green oval represents the end of one iteration and the red oval represents the end of one model run. The model has two internal counters: N t , which counts the number of agent time steps Δ t a , and N d , which counts the number of diffusion time steps Δ t d per Δ t a to be exactly AD . Time stepping The model updates its components using two different time steps. The first, Δ t d , governs the time evolution of the gradient grid. The second time step, Δ t a , governs updates to agent positions, orientations, and boundary strengths in the agent-based model (ABM) grid. This separation of time scales allows for efficient simulation while maintaining the accuracy of both the chemical diffusion dynamics and agent movement. The ratio of the two time steps is set as AD . Hence, during one iteration, we perform AD diffusion steps and one step of updating the agents and the ABM grid ( Figure 2 ). Calculate diffusion The model computes the time evolution of the chemical field s ( x, t ) on the gradient grid with a two-dimensional reaction-diffusion equation [16, 21, 37–40] with diffusion constant D s , decay rate d r and adsorption a d caused by the N agents at their ABM grid cells x n . As ABM grid cells are mapped to one gradient grid cell, the adsorption is normalized by that factor: This equation is discretized and numerically solved using the finite difference method FTCS (forward time-centered space) with Dirichlet boundary conditions, implemented using the Julia packages ParallelStencils and ImplicitGlobalGrid [ 41 ]. This method is commonly used for solving parabolic differential equations [ 42 – 44 ] and was specifically chosen because its time-explicit nature allows easier integration with the agent dynamics compared to implicit methods: The discretized equation 2 calculates the concentration s at each grid cell with indices ( i, j ) at time t + Δ t d , which changes from the value of through diffusion from or to neighboring grid cells, the decay of the chemical with rate d r , and the adsorption of the chemical with rate a d . The Dirac delta function ( δ ()) is only non-zero when x n is one of the ABM grid cells mapping to the gradient grid cell with position x i,j ( Fig. 1 ). The maximum time step size Δ t d to ensure numerical stability can be derived using von Neumann stability analysis [ 44 – 46 ] (see Section S.2 for more details): Agent update In every time step Δ t a , each agent updates its orientation and position. First, the orientation is updated according to: where the angle between two vectors is calculated as: and the normalization function N o is: The agents try to turn their current direction towards the negative gradient of the chemical field s ( x i ) at their position x i . The speed of this negative chemotactic alignment process is determined by the maximum gradient in the system through the normalization factor N o and a tunable turning rate parameter α , creating a relative response: For the case where α = 1, an agent at the location of the steepest gradient in the system will completely align its direction in one time step, while agents in regions with weaker gradients will align proportionally slower. When α < 1, the alignment speed is linearly scaled down by this factor. The term η i ( t ) models the stochastic component of the agent’s movement and alignment process and is drawn from a normal distribution with mean 0 ° and standard deviation η . Each agent’s position is updated according to: The agents move from their current position x i ( t ) in the direction θ i ( t + Δ t a ) with velocity v i , which is drawn from a normal distribution with mean v 0 and standard deviation v σ . Since x ( t + Δ t a ) is a vector of two continuous values, it is rounded to the nearest discrete grid cell on the ABM grid. Such an update scheme would be labeled as microscopic dry active matter [ 34 ], or in other terminology, modified active Brownian particles with auto-chemotactic alignment. Agent update depends on the time step Δ t a . Its value must strike a balance between two opposing constraints. If Δ t a is too small, the time discretization becomes overly fine, potentially introducing discretization artefacts: agents may move only one grid cell, or remain stationary, per step, which distorts movement dynamics [ 47 ]. Conversely, if Δ t a is too large, agents may traverse unrealistically long distances in a single step. This results in other discretization artefacts, such as unnaturally straight trajectories, which contradict experimental observations of cell motility [ 4 ]. In addition, rapid changes in the local environment, such as the number and identity of neighbouring agents or proximity to boundaries, impair biologically plausible agent-agent interactions. To balance these constraints, Δ t a was chosen to allow agents to move on average three ABM grid cells per time step in their current direction. This corresponds to a distance of 11.25 µ m, which is shorter than the typical cell length of a procyclic trypanosome (∼ 15–25 µ m [ 4 , 48 ]), thereby preserving local agent-agent interactions. For typical model parameters, the ratio ranges between 1 and 300. If necessary, the diffusion time step Δ t d is slightly reduced from its maximum allowable value to ensure that AD is always an integer. This guarantees a consistent number of diffusion time steps Δ t d per agent time step Δ t a throughout the simulation. Boundary Interactions The model simulates agent interactions with the colony boundary through reflective boundary conditions, whereby agents are reflected upon encountering a non-walkable grid cell on their path. Whilst real trypanosomes do not undergo physical reflection from colony boundaries, they typically remain at boundaries for only brief periods (a few seconds) before reorienting [ 4 ]. Therefore, reflective boundary conditions represent a justified simplification that captures the average behaviour at the agent numbers simulated here. The lattice constant Δ x and time step Δ t a are set to values so that agents move an average of three grid cells per time step. Since their movement paths can therefore interfere with non-walkable grid cells with k < 0, their movement paths are divided into segments of single grid cells. At each segment, the model checks for the presence of non-walkable grid cells with k < 0. If one is encountered, the agent is reflected. The reflection process approximates the local boundary morphology by constructing a circle with radius R t around the collision point (the first grid cell in an agent path segment with k < 0). The surface tangent is approximated using the two most distant non-boundary points neighbouring a boundary point. The agent is then reflected off this tangent such that the incoming angle β 1 equals the outgoing angle β 2 ( Fig. 3 ). Download figure Open in new tab Fig 3. Schematic representation of the boundary reflection mechanism. The boundary morphology is approximated at the collision point with a circle of radius R t , where its most distant surface points form a tangent. The agent is reflected with the outgoing angle β 2 equal to the incoming angle β 1 relative to the tangent. If an agent collides with a non-walkable grid cell on the colony boundary, the boundary weakens ( Fig. 4 ). This process is modelled through a change in the value k , which is initialized at all non-walkable grid cells with a negative number k 0 representing the grid strength. Upon collision, all grid cells with k < 0 within a circle of radius R c (distinct from R t ) centred on the colliding agent have their k values increased by 1. To address lattice anisotropy effects, not a perfect circle with radius R c was used, but rather a hybrid shape combining a circle and a diamond (see Section S.3 for more details). If k reaches zero, the non-walkable cells become walkable space and are removed from the boundary. Similar approaches have been successfully used to model bacterial colony growth [ 22 , 37 ]. Download figure Open in new tab Fig 4. Schematic representation of the boundary removal mechanism for an initial boundary strength k = − 3 and three subsequent agent collisions (a, b, c). Upon agent collision with a boundary, the boundary strength k increases by 1 in all boundary cells within radius R c . When k reaches zero after multiple collisions, boundary cells are converted to walkable space. Due to potentially complex local boundary morphology forming structures such as narrow channels, this collision-reflection-weakening process can occur up to 4 times in a single time step ( Fig. 2 ). Such complex morphology can also create cases where no unambiguous tangent can be determined; in these situations, a fallback mechanism is invoked where the agent is not reflected but simply reverses its original direction. To simulate a simple version of surface tension, a grid recovery rate g rc is implemented. This parameter decreases the absolute value of k in every boundary grid cell by a random integer drawn from a Poisson distribution with λ = g rc · Δ t a in each time step Δ t a until it reaches its initialization value k 0 , effectively strengthening the boundary over time if it is not further weakened by agent collisions. Cell Division Cell growth is modelled as exponential, based on population doubling times τ d observed in experiments [ 3 , 4 , 12 ]. From the doubling time, we derive the growth rate r as: In each time step Δ t a , each agent has a probability r · Δ t a of dividing and forming a new agent in a neighbouring empty grid cell. If no neighbouring empty grid cell is available, the agent cannot divide due to spatial constraints. The new agent is assigned a random orientation independent of the parent agent. Throughout the results and discussion sections, we primarily use the doubling time τ d rather than the growth rate r to describe growth dynamics, as it provides a more intuitive interpretation of population expansion. Model visualisation The visualisation integrates all three components of the model into a single composite image ( Figure 5 ). The agents are represented by a vector plot, where each arrow corresponds to one agent. These arrows indicate the current direction of movement, with additional colour-coding to distinguish agent orientations. Download figure Open in new tab Fig 5. Visualisation scheme of the model for default parameters ( Table 1 ) with L = x = y = 15 mm and after 9 hours: (a) The ∼ 5 × 10 5 agents are displayed as a vector plot where each agent is represented by an arrow pointing in its current direction. The direction is additionally colour-coded. Due to the number of agents, it is not possible to distinguish individual agents. (b) The gradient grid is visualized by a heat map where relative concentration values map to colours (colourbar on the right). (c) The ABM grid is displayed as a binary mask where non-colony regions ( k < 0) are shaded gray and regions within the colony ( k ≥ 0) are transparent. (d) The three components of the model are plotted together in one single plot containing all information. (e) Zoomed-in view of the complete plot to better visualize individual agents in their crowded environment. The agents are colour-coded according to their direction. The chemical gradient is displayed as a heat map with concentration values mapped to a colour scale. Since it is not yet clear which chemical is responsible for chemotactic alignment in trypanosome colonies and concentrations for the candidate chemicals have not been measured either, the absolute values of the concentration are less important than the relative concentration gradients in the system, which drive agent behaviour. Therefore, in all visualisations, the colourbar is auto-scaled to the maximum concentration present in each specific image. This approach allows us to visualize the relative concentration gradients consistently across all images using the same colour scaling, even though the absolute concentration values may differ substantially between different simulation conditions or time points. The ABM grid is visualized as a binary mask where non-colony regions appear gray while colony regions are transparent. These three visual elements are layered to provide a comprehensive view of the entire system, enabling intuitive observation of emerging patterns during simulation. Most images shown in the subsequent chapters do not display individual agents, as it is not possible to discern individual agents for simulations with realistic agent numbers, and most analysis in the following sections focuses primarily on the morphology of the ABM grid. Implementation The code for implementing the model and the simulations was designed with modularity and expandability as primary considerations. While this approach introduces some computational overhead, as some code parts are not optimized for maximum memory efficiency, it offers significant advantages for future development. The data structures have been chosen to allow a straightforward integration of additional model parameters or agent attributes in subsequent model iterations. The code is written in the programming language Julia [ 49 ], which aligns with the model’s design philosophy. Julia combines the readable, accessible syntax of high-level languages with performance comparable to traditional compiled languages like C and Fortran . This combination is particularly valuable for scientific computing applications where accessibility and runtime performance are critical. Julia’ s growing ecosystem of scientific libraries also provides valuable tools for efficient numerical computations, particularly for the reaction-diffusion equation. The code is available at https://github.com/scfischer/2026-Kuhn-et-al and will be archived at Zenodo upon publication. Simulation Modes The simulations can be run in two distinct modes: The first is the interactive mode, where only the current state of the simulation is stored in memory and the output is displayed with live animation using the Makie.jl plotting package [ 50 ]. In this mode, the parameter values can be dynamically adjusted during the simulation run to study their influence on model behaviour. The interactive mode can simulate up to 10 6 agents on modern CPUs (2025) in real-time, meaning one second in the model corresponds to one second or less in real time. The computation times of the model scale almost linearly with the number of agents and grid cells used. The second is the non-interactive mode, where the simulation runs with predetermined parameter values without graphical output. Here, the state of the simulation is saved at specified time points for later analyses. This mode is approximately twice as fast as the interactive mode and is primarily used for parameter sweeps in multiple parallel instances on workstations or clusters with sufficient memory per core. The saved data are subsequently used for analyses and visualisation. Parameter Space The model possesses a large parameter space, the parameters of which can be classified into four distinct categories: Parameters corresponding to measurable quantities that have been experimentally determined (e.g., growth rate r of trypanosomes in social motility assays). These experimentally measured values were adopted as default values in the model. Parameters corresponding to measurable quantities that have not yet been measured (e.g., adsorption rate a d , decay rate d r , and diffusion constant D s of the proposed chemotactic substance [ 11 , 13 ]). The reasons these remain unmeasured vary, but generally stem from either measurement difficulties or a lack of prior scientific motivation to quantify them. Parameters that do not directly correspond to measurable physical quantities (e.g., boundary strength k and grid recovery rate g rc ). These parameters arise from simplifications in the model. For example, the real colony boundary represents a complex hydrodynamic interface between a semi-fluid (agarose) and a fluid medium (nutrition solution containing trypanosomes), which can exchange liquid while remaining separated by surface tension. Although the properties of this interface could theoretically be measured, incorporating them would require modelling the interface with comparable complexity. The discretisation of time and space is also a model parameter; the values chosen/calculated for these are explained above and in Section S.1 . For parameters in categories (b) and (c), default values were determined through extensive testing, primarily in interactive mode. These values were selected to produce colony morphology dynamics similar to those observed and quantified in experiments [ 5 ]. A detailed analysis of their influence is shown in the results section. All parameters are shown in Table 1 together with their category and default values. There are two values given for the doubling time τ d as different experiments reported either doubling times of 9 h [ 3 , 12 ] or 20 h [ 4 ]. We tested both values. View this table: View inline View popup Download powerpoint Table 1. Model parameters categorised by type with their default values and units. The time steps Δ t a and Δ t d are dynamically calculated quantities that change depending on the diffusion coefficient D s , the mean agent velocity v 0 , the decay rate d r , and the lattice constant Δ x . Unless otherwise specified, all presented data are from simulated colonies starting with N = 2.5× 10 5 agents in a circular geometry with a diameter of 3 mm. This represents a similar trypanosome density to that observed in the well-quantified experiments from Kruger et al. [ 4 , 5 ], which serve as a blueprint for our model. In these experiments, approximately 10 6 agents are situated in colonies with diameters of 6 mm. Initial testing showed that this downscaling of the system by a factor of four did not change the overall behaviour but reduced computation times by a factor of approximately eight. This computational advantage aligns with other experimental observations, which demonstrated unchanged social motility activity in smaller colonies of N = 2×10 5 trypanosomes [ 12 ], hence its adoption in our simulations. Another performance-relevant quantity to determine before each simulation run is the edge length L of the simulated space. The computation time for the diffusion equation on the gradient grid scales with L 2 . Since we use Dirichlet boundary conditions, the colony cannot expand beyond the grid boundaries. Therefore, L should be large enough to prevent the colony from reaching the boundaries during the simulation period, yet as small as possible to maintain reasonable computation times. After testing, we determined that L = 15 mm for simulations running 10 h and L = 30 mm for simulations running 20 h satisfied these conditions well, and we used these values throughout our analyses. One important observation from our testing is that the variation between model runs with identical parameters is remarkably small. Given this high reproducibility and the significant computational demands, we opted to run each parameter set only once for the simulations presented below, thereby substantially reducing overall computation time without compromising the validity of our conclusions. Analysis The simulation results were analyzed by characterizing the ABM grid. We measured the area A of the colony at each time point as the number of colony grid points and normalised the increase in area with respect to the initial time point. Hence, for a given time point t ∗ and initial time point t 0 , we calculated the relative area increase of the simulated colony as . Furthermore, we counted the number of agents at each time point to obtain the cell number and calculated the sum of adsorbed material on the gradient grid. The cell density was calculated as . To quantify the colony morphology, we used surface roughness W , coefficient of variation CV , relative maximum peak height MaxP 1 and mean Fourier amplitude I k . These morphological metrics are based on an angular metric and a pair correlation metric and were previously introduced for trypanosome colonies [ 5 ]. Results Simulated colonies show two-phase expansion behaviour To establish baseline behaviour and validate our model against experimental observations, we simulated the system for 20 h using default parameter values ( Table 1 ) with two biologically relevant doubling times: τ d = 9 h ( Fig. 6 ) and τ d = 20 h ( Fig. 7 ), corresponding to different experimental conditions [ 3 , 4 , 12 ]. Download figure Open in new tab Fig 6. Time evolution of the colony morphology over 20 hours (2-hour intervals) with τ d = 9 h and default parameters ( Table 1 ). Initial noise in the gradient grid reflects random initialization. For details of visualisation see Materials and Methods. Download figure Open in new tab Fig 7. Time evolution of the colony morphology over 20 hours (2-hour intervals) with τ d = 20 h and default parameters ( Table 1 ). Initial noise in the gradient grid reflects random initialization. For details of visualisation see Materials and Methods. Both simulations exhibit a two-phase expansion: initial circular growth (0-2 h) followed by finger-like expansion. The faster-growing colony ( τ d = 9 h) develops approximately 50 fingers with more uniform morphology, while the slower case ( τ d = 20 h) produces approximately 45 fingers with greater variability in length. During the circular phase (0–2 h), the relative area increase is nearly identical across colonies, despite differences in doubling time ( Fig. 8 ). The number of cells increases exponentially throughout this phase. Decay and adsorption balance within the first hour. Afterwards, the amount of adsorbed material increases exponentially, mirroring cell growth. Cell density reaches a minimum during the first 10 h, coinciding with the morphological transition. All measures increase more rapidly in the faster-growing colony. Download figure Open in new tab Fig 8. Temporal evolution of colony properties in 30-minute intervals for two doubling times and default parameters ( Table 1 ). The results are for one simulation run each due to very low variability between runs (see Materials and Methods for details). (a) Relative area increase, (b) Cell number, (c) Adsorbed material, (d) Cell density. Note the density minimum at 4 h, coinciding with the morphological transition. We employed our established morphological metrics that are able to quantify the change from circular to finger-like growth by asserting values of zero to a colony growing perfectly circular and increasing to higher values for more finger-like growth ( Fig. 9 , [ 5 ]). Download figure Open in new tab Fig 9. Temporal evolution of morphological metrics in 30-minute intervals for two doubling times and default parameters ( Table 1 ). The results are for one simulation run each due to very low variability between runs (see Materials and Methods for details). (a) Surface roughness W , (b) Coefficient of variation CV , (c) Relative maximum peak height MaxP 1, and (d) Mean Fourier amplitude I k . All metrics clearly detect the initial circular growth phase, showing values at or near zero during the first two hours of the simulations. Subsequently, all metrics capture the transition from circular to finger-like growth as their values begin to increase. The coefficient of variation of the angular metric ( CV ) shows the highest sensitivity to this transition as it already increases at 1.5 h. The roughness ( W ) and the mean Fourier amplitude ( I k ) exhibit patterns of steady increase for both doubling times. In contrast, the CV and MaxP 1 metrics increase more slowly after 10 h. For τ d = 9 h, they eventually saturate. Both quantities particularly represent expansion that is perpendicular to the center point of the colony at t = 0 h [ 5 ]. After 10 h, this becomes less pronounced because single fingers begin to diverge from the initial straight paths and create a less uniform expansion front ( Fig. 6 ). In summary, the model exhibits the two-phase expansion behaviour observed in vitro . Our established metrics for in vitro Trypanosoma colonies successfully quantify the behaviour of the in silico colonies [ 5 ]. Simulations can reproduce experimental morphological metrics We compared our results with experimental data from [ 4 ] and [ 5 ]. Visual comparison shows similar morphologies but with more fingers in the simulated data compared to the experimental data. Additionally, our simulations show the perpendicular alignment of trypanosomes at the colony boundary that has been observed experimentally ( Fig 5 , [ 4 ]). We further conducted a quantitative comparison of the simulation results with respect to the experimental data for two different conditions [ 5 ]. In experiment 1, 10 6 cells were seeded. The colonies showed a relative increase in area of 1 after 20 h and developed 15 fingers. In experiment 2, twice the number of cells was seeded. The colonies showed a relative increase in area of 3.5 after 20 h and developed 30 fingers. Compared to this experimental data, the simulated colonies grow significantly faster, with a relative area increase of 7.5 for a doubling time of 20 h and 13 for a doubling time of 9 h, respectively ( Fig. 8a ). To enable a quantitative comparison of the morphological metrics despite the differences in relative area increase, we plotted the simulation results and the experimental data against relative area increase ( Fig. 10 and 11 ). The simulated colony with τ d = 20 h transitions to finger-like growth at a smaller relative area increase compared to the simulated colony with τ d = 9 h ( Fig. 10 ), indicating that increased growth rate not only accelerates colony expansion but also changes the nature of that expansion. For the experimental data, experiment 1 exhibits clearer fingering at a smaller relative area increase when compared to experiment 2 ( Fig. 11 ). We find that for experiment 2, the orders of magnitude of the simulation results agree well with the experimental data. For experiment 1, the morphological metrics are larger than for the simulations. Download figure Open in new tab Fig 10. Evolution of morphological metrics relative to area increase for default parameter values ( Table 1 ) with doubling times of τ d = 9 h and τ d = 20 h. The results are for one simulation run each due to very low variability between runs (see Materials and Methods for details). (a) Surface roughness W , (b) Coefficient of variation CV , (c) Relative maximum peak height MaxP 1, and (d) Mean Fourier amplitude I k . Download figure Open in new tab Fig 11. Evolution of morphological metrics relative to area increase for experimental data [ 47 ]. The dots indicate the mean and the error bars the standard deviation. (a) Surface roughness W , (b) Coefficient of variation CV , (c) Relative maximum peak height MaxP 1, and (d) Mean Fourier amplitude I k . Our results show an increase in colony growth speed and in number of fingers in simulations compared to experiments. Furthermore, we find a good agreement for the morphological metrics between the simulations and experiment 2. Parameter sensitivity analysis As the baseline behaviour was established, we performed an initial testing phase with the interactive simulation mode (see Materials and Methods for details). We then systematically varied model parameters that showed the most influence on the behaviour, or are mapped directly to open questions in the field. For all subsequent simulations, we used τ d = 9 h as the default value and ran each simulation for 10 h, as this duration was sufficient for the characteristic two expansion phases to emerge ( Fig. 6 ). We started our analysis with the agent’s movement, which is mostly determined by the orientation noise η and the turning rate α ( Equation (4) ). Orientation noise affects colony morphology The parameter η is the stochastic component of the direction alignment. The two edge cases η = 0 and η = 2 π represent perfect chemotactic alignment and complete random direction assignment at every time step, respectively. Neither extreme represents realistic microswimmer behaviour, as their active propulsion always has a stochastic components [ 34 ]. We varied η in 16 steps over a range of 0.0 – 2 π and analysed five parameter values in that range that represent the different model behaviours observed in that range ( Fig. 12 ). Download figure Open in new tab Fig 12. Colony morphology for five different values of the orientation noise η after 10 h with default parameters ( Table 1 ). For details of visualisation see Materials and Methods. For η = 0, the colony expands rapidly but irregularly, with small, irregular fingers of varying sizes. Increasing the noise slightly to η = 0.05 creates a more uniform gradient and more regular, finger-like expansion. At η = 0.2, the expansion pattern becomes much more regular, resembling the patterns observed in in vitro colonies. For higher noise values like η = 0.4, expansion is significantly slower, and fingers are less distinct. At η = 0.75, fingers disappear entirely, and expansion becomes exclusively circular and much slower. The area increases faster for lower noise values ( Fig. 13 ). Cell numbers and adsorbed material are nearly identical across all systems. The cell density increases for larger levels of noise. Download figure Open in new tab Fig 13. Temporal evolution of colony metrics over 10 h for default parameters ( Table 1 ) with five different noise values in 20-minutes intervals. The results are for one simulation run each due to very low variability between runs (see Materials and Methods for details). (a) Relative area, (b) Cell number, (c) Adsorbed material, (d) Cell density within the colony, (e) Surface roughness W , (f) Coefficient of variation CV , (g) Relative maximum peak height MaxP 1, and (h) Mean Fourier amplitude I k . Evolution of morphological metrics relative to area increase: (i) Surface roughness W , (j) Coefficient of variation CV , (k) Relative maximum peak height MaxP 1, and (l) Mean Fourier amplitude I k . The morphological metrics assume the highest values for orientation noise η = 0.2. W and I k show bigger differences when normalized for area increase, indicating a higher sensitivity to differences in relative area increase (panels i, l). In contrast, MaxP 1 only assumes values bigger than 0.1 for η = 0.2, which has very regular fingers perpendicular to the colony center. Hence, this value is highly sensitive in detecting perpendicular expansion behaviour, independent of area increase. When normalizing for relative area increase, the colony with η = 0.4 shows an early but less steep increase in all four metrics, indicating an early transition to finger growth but less pronounced fingers compared to the default parameters with η = 0.2. The CV shows constant values for η = 0.75, which is due to the slowly added, low number of grid cells to the colony, causing the grid anisotropy to become the dominant shaping factor. In summary, the orientation noise η has a huge impact on model behaviour, and colony development similar to in vitro colonies only occurs for a small range of parameter values. Turning rate affects expansion speed The turning rate α determines how quickly an agent responds to changes in the surrounding chemical concentration and aligns its orientation parallel to the negative gradient. The two edge cases represent instantaneous alignment ( α = 1.0) versus no alignment ( α = 0.0). We varied α in 16 steps over a range of 0.0–1.0 and analysed five parameter values in that range that represent the different model behaviours observed in that range ( Fig. 14 ). Download figure Open in new tab Fig 14. Colony morphology for five different values of the turning rate α after 10 h with default parameters ( Table 1 ). For details of visualisation, see Materials and Methods. The area increases more rapidly for higher turning rates, while cell number and adsorbed material remain nearly constant across all systems ( Fig. 15 ). The cell density decreases as α is increased. Analysis of the morphological metrics reveals a complex picture. The turning rate α = 0.1 shows the highest values for W and I k , slightly higher than α = 0.3. This occurs because both metrics are sensitive to finger formation and area increase, with finger formation being stronger for α = 0.1 but area increase being higher for α = 0.3. The CV and MaxP 1 metrics, which are less dependent on area increase and more sensitive to finger formation, show the highest values for α = 0.02, followed closely by the default value α = 0.1. Download figure Open in new tab Fig 15. Temporal evolution of colony metrics over 10 h for default parameters ( Table 1 ) with five different turning rates in 20-minutes intervals. The results are for one simulation run each due to very low variability between runs (see Materials and Methods for details). (a) Relative area increase, (b) Cell number, (c) Adsorbed material, (d) Cell density within the colony, (e) Surface roughness W , (f) Coefficient of variation CV , (g) Relative maximum peak height MaxP 1, and (h) Mean Fourier amplitude I k . Evolution of morphological metrics relative to area increase: (i) Surface roughness W , (j) Coefficient of variation CV , (k) Relative maximum peak height MaxP 1, and (l) Mean Fourier amplitude I k . When normalized for area increase, we observe that colonies with lower values of α grow fingers at a much lower area value than for α = 0.1 and especially for α = 0.3. However, for very low values such as α = 0.005, finger formation becomes so irregular that CV and MaxP 1 start to saturate after 5 h or a relative area increase of 1. We find that the duration for the circular expansion phase is independent of the value of α . The value of relative area increase for the onset of fingering, however, increases with increasing turning rate α . Small range of suitable parameter values for turning rate and orientation noise Since both the noise η and the turning rate α influence the movement behaviour of agents, we examined their coupled effects across 20 parameter combinations that represent the different model behaviours observed in that range ( Fig. 16 ). Download figure Open in new tab Fig 16. Colony morphology for five different values of the turning rate α (horizontal axis) and four different values of the noise η 0 (vertical axis) after 10 h with otherwise default parameters ( Table 1 ). For details of visualisation see Materials and Methods. Within the chosen parameter range of η from 0.1 to 0.4 and α from 0.005 to 0.3, the turning rate α has a greater impact on colony morphology. The noise η exhibits an important effect: if it is too high, lattice anisotropy becomes a defining factor in the emerging morphology. However, η cannot be too low either, as this leads to more irregular finger formation, especially pronounced for low values of α . From visual inspection, the parameter space where expansion behaviour is similar to experiments is quite small at η ≈ 0.2 and α ≈ 0.02–0.1, a range we have already analysed ( Fig. 15 ). Grid and boundary properties affect expansion speed After analysing the movement properties of the agents, we focused on the boundary properties in order to understand how such environmental factors impact the colony expansion. We varied the grid strength k 0 from −1000 to −25000 and the grid recovery rate g rc from 0.0 to 25.0 and analysed five parameter values each in those ranges that represent the different model behaviours observed ( Fig. 17 ). To make things more intuitive to understand, we always use the absolute value of | k 0 | in the following sections to avoid confusion regarding the negative values for k 0 . Download figure Open in new tab Fig 17. First row: Colony morphology for five different values of the grid strength | k 0 | after 10 h with default parameters ( Table 1 ) Second row: Colony morphology for five different values of the grid recovery rate g rc after 10 h with default parameters ( Table 1 ). For details of visualisation see Materials and Methods. Both a lower absolute grid strength | k 0 | and a lower grid recovery rate g rc cause faster expansion. However, whilst a lower g rc also causes the expansion behaviour to change to circular growth, a lower | k 0 | does not appear to change the expansion pattern. Similarly, low values of g rc of 0.0 and 2.0 show much higher relative area increase ( Fig. 18 ), whereas all other changes in both parameters only steadily increase or decrease the relative area. Cell number and absorbed material are not shown as the behaviour is the same as for the default parameters. Cell density reflects the change in relative area. Download figure Open in new tab Fig 18. Temporal evolution of relative area increase over 10 h for default parameters ( Table 1 ) in 20-minute intervals with (a) five different grid recovery rate values g rc and (b) five different absolute grid strength values | k 0 |. The results are for one simulation run each due to very low variability between runs (see Materials and Methods for details). All morphological metrics relative to time for different values of | k 0 | show equidistant spacing with a later increase in values but similar slopes ( Fig. 19 ). When normalized by relative area increase, the trend becomes clearer. Except at very high grid strength values of | k 0 |, all morphological metrics exhibit a marked increase within the same relative area interval of three to four. Download figure Open in new tab Fig 19. Evolution of morphological metrics over 10 h for default parameters ( Table 1 ) in 20-minutes intervals for different grid strength values | k 0 | relative to time (first row) and relative to area increase (second row); for different grid recovery rate values g rc relative to time (third row) and relative to area increase (fourth row). The results are for one simulation run each due to very low variability between runs (see Materials and Methods for details). Something very similar can be observed for the grid recovery rate g rc , where the time point for the onset of fingering increases with increasing g rc . For the evolution with respect to relative area increase, all curves align except for the very low values of g rc = 0.0 and g rc = 2.0. We conclude that there is a wide parameter corridor for k 0 and g rc that changes the temporal onset of fingering, i.e. the expansion speed but not the nature of the expansion. Diffusion coefficient affects colony morphology and expansion speed Following the previous results demonstrating that boundary properties have less influence on colony expansion behaviour than the movement properties of agents, we turned our attention to another major factor that influences agent movement: the diffusion of the chemotactically active substance. Various substances (pH/protons, glucose, cAMP, exosomes) [ 10 – 13 ] have been proposed as potential mediators of a chemotactic alignment in trypanosomes [ 11 ]. However, it remains unclear which specific substance could mediate the chemotactic alignment, as all are present in the experimental system. The diffusion coefficients of these candidate substances at room temperature in water are well established in the literature, ranging from 9× 10 −9 m 2 /s (pH/ H + ), 6.5× 10 −10 m 2 /s (glucose) [ 51 ], 4× 10 −10 m 2 /s (cAMP) [ 52 , 53 ] to 7× 10 −12 –1.5× 10 −12 m 2 /s (exosomes, depending on size) [ 10 , 54 – 56 ]. Although the experimental setup of trypanosome colonies involves a complex liquid medium that consists of multiple components [ 4 ] placed on top of an agarose gel, previous studies have shown that the diffusion coefficients of various materials in such media differ only slightly from those in pure water [ 57 , 58 ]. Therefore, we can reasonably use the diffusion coefficients measured in water as reference points. We varied the diffusion coefficient D s in our simulations across the range of biologically relevant values and studied its influence on system behaviour ( Fig. 20 ). Download figure Open in new tab Fig 20. Colony morphology for ten different values of the diffusion coefficient D s after 10 h with default parameters ( Table 1 ). For details of visualisation see Materials and Methods. The diffusion coefficient has a tremendous influence on the expansion behaviour. For the lowest value of D s = 1× 10 −12 m 2 /s, diffusion is so slow that no colony-wide concentration gradient decreasing from the center to the boundaries is established within 10 h. Instead the highest concentration is found at the boundaries of the colony. The resulting colony morphology is anisotropic and does not exhibit a finger-like pattern. For higher diffusion coefficients, a colony-wide gradient can be established, with the highest concentration in the central part of the colony decreasing outward. This gradient becomes more and more pronounced the higher the diffusion coefficient D s becomes. The resulting morphologies change steadily toward the very regular finger-like patterns observed at D s = 7× 10 −11 m 2 /s. For higher values up to 8.7× 10 −10 m 2 /s, expansion becomes even faster, but fingers start to branch more and also merge, forming a more circular and less fractal expansion front. For even higher values of 5× 10 −9 m 2 /s and beyond, expansion becomes slower again. Thick protrusions emerge that are not comparable to the finger-like patterns observed in experiments. Starting at D s = 3.7× 10 −10 m 2 /s, the adsorbed chemical diffuses significantly beyond the colony fingers, which can be seen in the changing colours outside of the colony boundaries ( Fig. 20 ), a trend that increases for larger D s values. At D s = 9× 10 −9 m 2 /s, the chemical and its decreasing concentration gradient reach far beyond the colony boundaries. The colony behaviour is the same for the first two hours ( Fig. 21 ). After that, the area increases faster for higher diffusion coefficients until D s = 8.7 × 10 −10 m 2 /s. Then, the area increases more slowly for higher values of D s . The cell number and adsorbed material are similar for all values of D s until D s = 7.0 × 10 −11 m 2 /s and then decrease with increasing D s . The cell density mirrors the relative area increase. Download figure Open in new tab Fig 21. Temporal evolution of colony metrics over 10 h for default parameters ( Table 1 ) with five representative diffusion coefficient D s values in 20-minutes intervals. The results are for one simulation run each due to very low variability between runs (see Materials and Methods for details). (a) Relative area increase, (b) Cell number, (c) Total amount of material in the gradient grid, (d) Cell density within the colony, (e) Surface roughness W , (f) Coefficient of variation CV , (g) Relative maximum peak height MaxP 1, and (h) Mean Fourier amplitude I k . Evolution of morphological metrics relative to area increase: (i) Surface roughness W , (j) Coefficient of variation CV , (k) Relative maximum peak height MaxP 1, and (l) Mean Fourier amplitude I k . Details of visualisation in Sec.. Analysis of the morphological metrics reveals a complex picture. The default diffusion coefficient D s = 7.0× 10 −11 m 2 /s shows the highest values for CV and MaxP 1, slightly higher than D s = 2.1× 10 −11 m 2 /s. This occurs because both metrics are sensitive to finger formation independent of area increase, with finger formation being slightly stronger for D s = 2.1× 10 −11 m 2 /s. This can be seen in W and I k , which are both sensitive to finger formation and area increase, showing the highest values for D s = 8.7× 10 −10 m 2 /s, slightly higher than for D s = 2.1× 10 −11 m 2 /s, as finger formation is stronger in the latter, while area increase is stronger in the former. When normalised for area increase, all metrics behave very similarly. We observe that colonies with lower values of D s start forming fingers at a lower relative area. However, for the low value of D s = 6.0× 10 −12 m 2 /s, all morphological metrics also indicate a beginning saturation in finger formation after a relative area increase of approximately 3. The system with default diffusion coefficient D s = 7.0× 10 −11 m 2 /s shows a later onset of non-circular expansion but a similar increase in the metrics, indicating similar rapid expansion in finger-like patterns. For higher values of D s , all metrics show a later onset of increase and a slower slope, indicating a less finger-dominated expansion pattern. In summary, the diffusion coefficient has a clear effect on the onset of fingering with respect to relative area increase and on finger morphology. Our parameter sensitivity analysis shows that the different parameters analysed exhibit distinct effects on model behaviour. Movement parameters η and α demonstrate a narrow parameter window for producing behaviour similar to experiments. Boundary property parameters k 0 and g rc primarily affect expansion speed with minimal influence on expansion pattern. The diffusion coefficient D s exerts a substantial influence on expansion dynamics. Low values result in slow, irregular colony expansion. Intermediate values enhance both expansion speed and regularity. For high values, this trend reverses such that the colonies expand more slowly and regularly but exhibit progressively reduced finger formation. Discussion Our agent-based model demonstrates that complex colony-level patterns in trypanosome social motility can emerge from the combination of single cell motility and proliferation, interactions of the cells with the colony boundaries and negative autochemotaxis. Our simulation results are sensitive to the parameters for growth dynamics and single cell motility. The two simulated doubling times demonstrate a non-linear relationship between cell growth and area increase. The underlying mechanism remains unclear, though relative area growth may be primarily driven by agents at the colony surface, where density appears similar across both systems. This requires further investigation through comparison with experimental data using cell lines with different doubling times [ 3 , 4 , 12 ] and detailed analysis of agent position and density distributions over time within the model. The orientation noise η and turning rate α demonstrate a narrow parameter window where finger formation is most pronounced, regular, and most similar to experiments. Such behaviour represents a remarkable result. It suggests that finger formation emerges as a direct consequence of the specific properties of the directional random walk of individual agents. Our results coincide with previous experimental findings where modifying flagellum activity and motility properties of individual trypanosomes can suppress social motility [ 1 , 3 ]. Although extensive work has examined individual trypanosome movement properties both experimentally [ 29 ] and through modelling [ 27 , 28 ], direct comparison with our findings remains challenging due to fundamental differences between three-dimensional free environments and our crowded two-dimensional system that constrains movement. It has been demonstrated that agent trajectories can be measured in the in vitro system [ 4 , 11 ]. Such data is very suitable for a future direct comparison with agent trajectories extracted from the model. The boundary parameters grid strength k 0 and grid recovery rate g rc exhibit a wide “parameter corridor” where changes affect only expansion speed without altering the fundamental nature of colony morphologies. This robustness suggests that the specific mechanical properties of the colony-agarose interface may be less critical than previously thought [ 4 ], with qualitative expansion behaviour being primarily determined by the interplay between chemotaxis and movement dynamics. However, our current boundary implementation is quite simplistic and could be expanded to much more complex hydrodynamic models that have been implemented for bacterial systems like Pseudomonas [ 23 ], albeit at much smaller scales or in a continuous manner where the colony boundary is treated as a free boundary problem [ 21 , 40 ]. Significant computational optimisation and modification would be needed to make hydrodynamic boundary models tractable and to couple them to an agent-based model with the number of agents we have used here. Compared to our phenomenological implementation of boundaries, such models also require measurements that have not yet been done to quantify the defining physical parameter values of the hydrodynamic interface, and all of them assume some sort of lubrication caused by the microswimmers [ 24 , 59 ] themselves that changes the local properties of the boundary. Whether trypanosomes secrete substances or lubricate their environment in a significant manner is currently only speculative. The diffusion coefficient D s emerges as the most influential parameter in our model, dramatically affecting both the speed and pattern of colony expansion. For low values of D s = 1.0–6× 10 −12 m 2 /s, finger formation is absent and colony growth is substantially reduced, likely because information cannot spread sufficiently rapidly across the colony to establish a global chemical gradient pointing outward toward the boundaries. Conversely, very high diffusion coefficients lead to more circular rather than finger-like growth. At these high values, signals diffuse rapidly and concentration gradients become shallow but far-reaching, extending beyond the colony boundaries. Such conditions inhibit finger formation because the highest gradient is no longer localised at the finger tips, where it guides directional growth, but instead decreases monotonically with distance from the colony centre, resulting in a unified expansion front without finger formation. The range of simulated diffusion coefficients (10 −12 to 10 −9 m 2 /s) spans the known values for several proposed signalling agents, including pH (H + ), cAMP, glucose, and small exosomes [ 3 , 10 , 12 ]. However, the model’s behaviour most similar to experiments occurs within the range of D s = 2× 10 −11 to 10 −10 m 2 /s. Using the Einstein-Stokes equation: where k B is the Boltzmann constant, T = 293.15 K is the absolute temperature, and η = 0.001 kg/(m·s) is the dynamic viscosity of water at T , the corresponding hydrodynamic radius r of the molecule can be estimated to be between 2–10 nm ( Fig. 22 ). Download figure Open in new tab Fig 22. Diffusion coefficient as a function of hydrodynamic radius of diffusion particles in water according to the Einstein-Stokes equation. Model behaviour similar to experiments occurs at diffusion coefficients between 2× 10 −11 and 10 −10 m 2 /s, corresponding to hydrodynamic radii between 2–10 nm. Using established conversion estimations for molecular weights [ 60 ], such radii correspond to small proteins with molecular weights between 12.1–1690 kDa. Therefore, our results suggest searching for signalling proteins/macromolecules/lipids in that weight range. A quantitative comparison between our simulation results and experimental data shows that in the simulations, the colony growth speed and the number of fingers are larger. In contrast, the morphological metrics agree well between simulations and experiments. Hence, while the overall growth rate differs, the pattern formation dynamics remain comparable. Our parameter sensitivity analysis further indicates that adjusting the boundary strength and recovery parameters could improve the quantitative agreement. Previous studies have shown that social motility can also occur with chiral fingers — a behaviour that is currently not captured by our model. Similarly, the inhibition of finger formation, for example, by a neighbouring E. coli colony, is not yet represented. Addressing these phenomena would require formulating hypotheses about the underlying mechanisms that could alter local motility, such as asymmetric propulsion or interspecies signalling, and subsequently incorporating these processes into the model. Conclusion Our model reproduces the patterning characteristics of social motility of single colonies and demonstrates that individual movement properties of trypanosomes are critical for the emergence of finger-like patterns. Our main prediction is that boundary interactions have to be coupled with negative auto-chemotaxis and that previously proposed auto-chemotactically active substances cannot be directly responsible for social motility and likely represent correlation rather than causation. To advance this research, we have provided several testable hypotheses for experimental validation, including the motility statistics of trypanosomes, the size of potential signalling agents, and the robustness with respect to mechanical boundary properties. Supporting information S.1 Area of Cells The projected two-dimensional area of trypanosomes during a social motility assay was determined using scanning electron microscope images from [ 4 ]. The cell body was manually segmented with Fiji, and the pixel count was used to calculate an approximate area of 18.0 µ m 2 . Download figure Open in new tab Fig S.1. Scanning electron micrograph of a single trypanosome fixed while swimming on the surface of an agarose gel. The manually determined outline of the cell body is shown in orange. A and P indicate the anterior and posterior ends, respectively. Scale bar: 5 µ m. Image adapted from [ 4 ]. S.2 Von Neumann Stability Analysis The von Neumann stability analysis for standard diffusion equations solved with finite difference schemes has been extensively documented in the literature [ 44 – 46 ]. Here, we focus specifically on the stability implications when extending the analysis to our modified diffusion equation that includes additional decay and adsorption terms: where D s is the diffusion coefficient, d r represents the decay rate, s f is the scaling factor between the ABM grid and gradient grid, a d is the adsorption rate caused by the N agents at their ABM grid cells x n , δ is the Dirac delta function. For stability analysis, we consider the worst-case scenario where all ABM grid points are occupied by agents. In this case, the adsorption term becomes a constant at every gradient grid point. The equation can be written in the form: Although this equation contains a non-homogeneous term C , the stability of the numerical scheme depends only on the homogeneous part of the equation. This can be seen by examining the round-off error , which is the difference between the numerical solution and the solution with infinite machine precision: Both solutions must satisfy our discretized equation: If we let L ∗ represent our discretized operator, then for : And similarly for the numerical solution: By the linearity of the operator L ∗ , the round-off error satisfies: Therefore, the round-off error satisfies the homogeneous equation, allowing us to apply the standard von Neumann stability analysis for the explicit forward time-centered space scheme. The derivation follows the same steps as in [ 45 ], but extended to 2D (see [ 61 ] for the general 2D case) and includes an additional decay term. Therefore, we can write as a Fourier series: We examine the amplification factor which represents the relative change of error between iterations: Using the identity e ik Δ x + e − ik Δ x = 2 cos( k Δ x ), this can be rearranged to: Since cos( θ ) − 1 ≤ 0 for all θ , the most restrictive condition occurs when cos( k x Δ x ) = cos( k y Δ y ) = −1 [ 61 ], giving: For stability, we require | G | ≤ 1, which leads to: The right inequality is always satisfied, and from the left inequality we get: Which simplifies to: As in our simulation setup Δ x = Δ y , this can be further simplified and rearranged to our final stability condition on Δ t : For the standard diffusion equation (when d r = 0), this reduces to the familiar condition [ 61 ]. Importantly, the presence of the decay term d further restricts the maximum allowable time step, requiring an even smaller Δ t to maintain numerical stability in our simulations. S.3 Lattice Anisotropy On-lattice cluster growth models can suffer from a phenomenon where the lattice symmetry becomes imprinted on the emerging structure, even with a completely random selection of growth sites [ 62 – 65 ]. In square lattices, this manifests as preferential growth along the major lattice axes (x and y) and their diagonals compared to other directions. This undesired anisotropy varies depending on the specific model and mitigation strategies employed. Most solutions require more complex surface site selection algorithms during the growth process to avoid anisotropic growth [ 63 , 64 , 66 , 67 ]. Our on-lattice model of trypanosome colony geometry can be understood in a broader sense as a cluster growth model, although the colony does not “grow” in the traditional sense. Instead, swimming trypanosomes push the colony boundary outward, with each agent influencing multiple grid cells inside the circular collision area when colliding with the boundary. Nevertheless, cell growth remains a significant factor enabling colony expansion. To investigate whether colony anisotropies do emerge in our model and how they might be mitigated, we conducted a simulation with randomly moving agents and high boundary strength (similar to the noise reduction parameter m in Eden-like models, which enhance growth anisotropies [ 62 ]). The simulations were normalized to produce a similar relative area increase after 10 hours. Anisotropy emerged when using a circular collision between agents and the boundary ( Fig. S.2 ). The colony preferentially grows along the x and y lattice axes, visible as four distinct peaks in the angular metric. Changing the collision geometry to a diamond shape reversed this trend completely, causing the colony to grow preferentially along the diagonals of the lattice. Conversely, a square-shaped collision area dramatically increased the initially observed preferential growth along the x and y axes. Based on these findings, we implemented a hybrid collision area combining both circular and diamond geometries. This approach effectively balanced the competing anisotropies, resulting in more isotropic colony growth that better matches experimental observations. Download figure Open in new tab Fig S.2. Comparison of different collision area geometries and their effect on colony growth anisotropy, normalized for similar relative area gain over 10 hours. The angular metric shows the distinct periodic anisotropies in the surface morphology. (a) Using a circular collision area results in preferential growth along the x and y axes. (b) A diamond-shaped collision area reverses this trend, favouring diagonal growth. (c) A square-shaped collision area further enhances growth along the x and y axes. (d) Hybrid collision geometry combines circle and diamond shapes to produce the most isotropic growth. Acknowledgements This project was funded by the Priority programme 2332 ‘Physics of Parasitism’ of the German Research Foundation (DFG). References 1. ↵ Oberholzer M , Lopez MA , McLelland BT , Hill KL . Social Motility in African Trypanosomes . PLoS Pathogens . 2010 ; 6 ( 1 ): e1000739 . doi: 10.1371/journal.ppat.1000739 . OpenUrl CrossRef PubMed 2. ↵ Imhof S , Knüsel S , Gunasekera K , Vu XL , Roditi I. Social Motility of African Trypanosomes Is a Property of a Distinct Life-Cycle Stage That Occurs Early in Tsetse Fly Transmission . PLOS Pathogens . 2014 ; 10 ( 10 ): e1004493 . doi: 10.1371/journal.ppat.1004493 . OpenUrl CrossRef PubMed 3. ↵ Oberholzer M , Saada EA , Hill KL . Cyclic AMP Regulates Social Behavior in African Trypanosomes . mBio . 2015 ; 6 ( 3 ): e01954 – 14 . doi: 10.1128/mBio.01954-14 . OpenUrl CrossRef PubMed 4. ↵ Krüger T , Maus K , Kreß V , Meyer-Natus E , Engstler M. Single-cell motile behaviour of Trypanosoma brucei in thin-layered fluid collectives . The European Physical Journal E . 2021 ; 44 ( 3 ): 37 . doi: 10.1140/epje/s10189-021-00052-7 . OpenUrl CrossRef 5. ↵ Kuhn A , Krüger T , Schüttler M , Engstler M , Fischer SC . Quantification of Trypanosoma brucei social motility indicates different colony growth phases . Journal of The Royal Society Interface . 2024 ; 21 ( 221 ): 20240469 . doi: 10.1098/rsif.2024.0469 . OpenUrl CrossRef PubMed 6. ↵ Langousis G , Hill K. Motility and more: the flagellum of Trypanosoma brucei . Nature Reviews Microbiology . 2014 ; 12 : 505 – 518 . doi: 10.1038/nrmicro3274 . OpenUrl CrossRef PubMed 7. Bargul J , Jung J , Mcodimba F , Omogo CO , Adung’a V , Krüger T , et al. Species-Specific Adaptations of Trypanosome Morphology and Motility to the Mammalian Host . PLoS Pathogens . 2016 ; 12 . doi: 10.1371/journal.ppat.1005448 . OpenUrl CrossRef PubMed 8. Sun SY , Kaelber J , Chen M , Dong X , Nematbakhsh Y , Shi J , et al. Flagellum couples cell shape to motility in Trypanosoma brucei . Proceedings of the National Academy of Sciences of the United States of America . 2018 ; 115 : E5916 – E5925 . doi: 10.1073/pnas.1722618115 . OpenUrl Abstract / FREE Full Text 9. ↵ Hill K. Parasites in motion: flagellum-driven cell motility in African trypanosomes . Current opinion in microbiology . 2010 ; 13 4: 459 – 65 . doi: 10.1016/j.mib.2010.05.015 . OpenUrl CrossRef PubMed 10. ↵ Eliaz D , Kannan S , Shaked H , Arvatz G , Tkacz ID , Binder L , et al. Exosome secretion affects social motility in Trypanosoma brucei . PLOS Pathogens . 2017 ; 13 ( 3 ): e1006245 . doi: 10.1371/journal.ppat.1006245 . OpenUrl CrossRef PubMed 11. ↵ DeMarco SF , Saada EA , Lopez MA , Hill KL . Identification of Positive Chemotaxis in the Protozoan Pathogen Trypanosoma brucei . mSphere . 2020 ; 5 ( 4 ): 10.1128/msphere.00685–20 . doi: 10.1128/msphere.00685-20 . OpenUrl CrossRef 12. ↵ Shaw S , Knüsel S , Abbühl D , Naguleswaran A , Etzensperger R , Benninger M , et al. Cyclic AMP signalling and glucose metabolism mediate pH taxis by African trypanosomes . Nature Communications . 2022 ; 13 ( 1 ): 603 . doi: 10.1038/s41467-022-28293-w . OpenUrl CrossRef PubMed 13. ↵ Shaw S , Roditi I. The sweet and sour sides of trypanosome social motility . Trends in Parasitology . 2023 ; 39 ( 4 ): 242 – 250 . doi: 10.1016/j.pt.2023.01.001 . OpenUrl CrossRef PubMed 14. ↵ Kearns DB , Losick R. Swarming motility in undomesticated Bacillus subtilis . Molecular Microbiology . 2003 ; 49 ( 3 ): 581 – 590 . doi: 10.1046/j.1365-2958.2003.03584.x . OpenUrl CrossRef PubMed Web of Science 15. Swiecicki JM , Sliusarenko O , Weibel DB . From swimming to swarming: Escherichia coli cell motility in two-dimensions . Integrative biology : quantitative biosciences from nano to macro . 2013 ; 5 ( 12 ): 1490 – 1494 . doi: 10.1039/c3ib40130h . OpenUrl CrossRef PubMed 16. ↵ Kawasaki K , Mochizuki A , Matsushita M , Umeda T , Shigesada N. Modeling Spatio-Temporal Patterns Generated byBacillus subtilis . Journal of Theoretical Biology . 1997 ; 188 ( 2 ): 177 – 185 . doi: 10.1006/jtbi.1997.0462 . OpenUrl CrossRef PubMed Web of Science 17. ↵ Bees MA , Andresén P , Mosekilde E , Givskov M. Quantitative Effects of Medium Hardness and Nutrient Availability on the Swarming Motility of Serratia liquefaciens . Bulletin of Mathematical Biology . 2002 ; 64 ( 3 ): 565 – 587 . doi: 10.1006/bulm.2002.0287 . OpenUrl CrossRef PubMed 18. ↵ Wu Y , Jiang Y , Kaiser D , Alber M. Social Interactions in Myxobacterial Swarming . PLOS Computational Biology . 2007 ; 3 ( 12 ): e253 . doi: 10.1371/journal.pcbi.0030253 . OpenUrl CrossRef PubMed 19. F Copeland M, B Weibel D. Bacterial swarming: a model system for studying dynamic self-assembly . Soft Matter . 2009 ; 5 ( 6 ): 1174 – 1187 . doi: 10.1039/B812146J . OpenUrl CrossRef PubMed Web of Science 20. ↵ Bonachela JA , Nadell CD , Xavier JB , Levin SA . Universality in Bacterial Colonies . Journal of Statistical Physics . 2011 ; 144 ( 2 ): 303 – 315 . doi: 10.1007/s10955-011-0179-x . OpenUrl CrossRef 21. ↵ Kozlovsky Y , Cohen I , Golding I , Ben-Jacob E. Lubricating bacteria model for branching growth of bacterial colonies . Physical Review E . 1999 ; 59 ( 6 ): 7025 – 7035 . doi: 10.1103/PhysRevE.59.7025 . OpenUrl CrossRef PubMed 22. ↵ Ben-Jacob E , Levine H. Self-engineering capabilities of bacteria . Journal of The Royal Society Interface . 2006 ; 3 ( 6 ): 197 – 214 . doi: 10.1098/rsif.2005.0089 . OpenUrl CrossRef PubMed Web of Science 23. ↵ Du H , Xu Z , Shrout JD , Alber M. Multiscale modeling of Pseudomonas Aeruginosa swarming . Mathematical Models and Methods in Applied Sciences . 2011 ; 21 (supp01): 939 – 954 . doi: 10.1142/S0218202511005428 . OpenUrl CrossRef PubMed 24. ↵ Bru JL , Kasallis SJ , Zhuo Q , Høyland-Kroghsbo NM , Siryaporn A. Swarming of P. aeruginosa: Through the lens of biophysics . Biophysics Reviews . 2023 ; 4 ( 3 ): 031305 . doi: 10.1063/5.0128140 . OpenUrl CrossRef PubMed 25. ↵ Tian M , Wu Z , Zhang R , Yuan J. A new mode of swimming in singly flagellated Pseudomonas aeruginosa . Proceedings of the National Academy of Sciences . 2022 ; 119 ( 14 ): e2120508119 . doi: 10.1073/pnas.2120508119 . OpenUrl CrossRef PubMed 26. ↵ Qian C , Wong CC , Swarup S , Chiam KH . Bacterial tethering analysis reveals a “run-reverse-turn” mechanism for Pseudomonas species motility . Applied and Environmental Microbiology . 2013 ; 79 ( 15 ): 4734 – 4743 . doi: 10.1128/AEM.01027-13 . OpenUrl Abstract / FREE Full Text 27. ↵ Alizadehrad D , Krüger T , Engstler M , Stark H. Simulating the Complex Cell Design of Trypanosoma brucei and Its Motility . PLOS Computational Biology . 2015 ; 11 ( 1 ): e1003967 . doi: 10.1371/journal.pcbi.1003967 . OpenUrl CrossRef PubMed 28. ↵ Overberg FA , Jamshidi Khameneh N , Krüger T , Engstler M , Gompper G , Fedosov DA . Modelling motility of Trypanosoma brucei . PLOS Computational Biology . 2025 ; 21 ( 5 ): e1013111 . OpenUrl 29. ↵ Wheeler RJ . Use of chiral cell shape to ensure highly directional swimming in trypanosomes . PLOS Computational Biology . 2017 ; 13 ( 1 ): e1005353 . doi: 10.1371/journal.pcbi.1005353 . OpenUrl CrossRef PubMed 30. ↵ Tropea C , Bleckmann H , editors Babu SB , Schmeltzer C , Stark H. Swimming at Low Reynolds Number: From Sheets to the African Trypanosome . In: Tropea C , Bleckmann H , editors. Nature-Inspired Fluid Mechanics: Results of the DFG Priority Programme 1207 “Nature-inspired Fluid Mechanics” 2006-2012. Berlin, Heidelberg : Springer ; 2012 . p. 25 – 41 . 31. Elgeti J , Winkler RG , Gompper G. Physics of Microswimmers - Single Particle Motion and Collective Behavior . Reports on Progress in Physics . 2015 ; 78 ( 5 ): 056601 . doi: 10.1088/0034-4885/78/5/056601 . OpenUrl CrossRef PubMed 32. Theers M , Westphal E , Qi K , G Winkler R , Gompper G. Clustering of microswimmers: interplay of shape and hydrodynamics . Soft Matter . 2018 ; 14 ( 42 ): 8590 – 8603 . doi: 10.1039/C8SM01390J . OpenUrl CrossRef PubMed 33. ↵ Bárdfalvy D , Škultéty V , Nardini C , Morozov A , Stenhammar J. Collective motion in a sheet of microswimmers . Communications Physics . 2024 ; 7 ( 1 ): 1 – 8 . doi: 10.1038/s42005-024-01587-9 . OpenUrl CrossRef 34. ↵ Shaebani MR , Wysocki A , Winkler RG , Gompper G , Rieger H. Computational models for active matter . Nature Reviews Physics . 2020 ; 2 ( 4 ): 181 – 199 . OpenUrl 35. ↵ Schaar K , Zöttl A , Stark H. Detention Times of Microswimmers Close to Surfaces: Influence of Hydrodynamic Interactions and Noise . Physical Review Letters . 2015 ; 115 ( 3 ): 038101 . doi: 10.1103/PhysRevLett.115.038101 . OpenUrl CrossRef PubMed 36. ↵ Elgeti J , Gompper G. Microswimmers near surfaces . The European Physical Journal Special Topics . 2016 ; 225 ( 11 ): 2333 – 2352 . doi: 10.1140/epjst/e2016-60070-6 . OpenUrl CrossRef 37. ↵ Ben-Jacob E , Schochet O , Tenenbaum A , Cohen I , Czirók A , Vicsek T. Generic modelling of cooperative growth patterns in bacterial colonies . Nature . 1994 ; 368 ( 6466 ): 46 – 49 . doi: 10.1038/368046a0 . OpenUrl CrossRef PubMed Web of Science 38. J Marsden E , Valeriani C , Sullivan I , E Cates M , Marenduzzo D. Chemotactic clusters in confined run-and-tumble bacteria: a numerical investigation . Soft Matter . 2014 ; 10 ( 1 ): 157 – 165 . doi: 10.1039/C3SM52358F . OpenUrl CrossRef PubMed 39. Blanchard AE , Lu T. Bacterial social interactions drive the emergence of differential spatial colony structures . BMC Systems Biology . 2015 ; 9 ( 1 ): 59 . doi: 10.1186/s12918-015-0188-5 . OpenUrl CrossRef PubMed 40. ↵ Giverso C , Verani M , Ciarletta P. Emerging morphologies in round bacterial colonies: comparing volumetric versus chemotactic expansion . Biomechanics and Modeling in Mechanobiology . 2016 ; 15 ( 3 ): 643 – 661 . doi: 10.1007/s10237-015-0714-9 . OpenUrl CrossRef PubMed 41. ↵ Omlin S , Räss L , Utkin I. Distributed Parallelization of xPU Stencil Computations in Julia ; 2022 . 42. ↵ Roache PJ . Computational fluid dynamics . Computational Fluid Dynamics . 1976 ;. 43. Anderson DA , Tannehill JC , Pletcher RH . Computational fluid mechanics and heat transfer. BOOK . Taylor & Francis ; 1997 . 44. ↵ Blazek J. Computational Fluid Dynamics: Principles and Applications . Butterworth-Heinemann ; 2015 . 45. ↵ John D , Anderson J. In: Computational fluid dynamics: the basics with applications . McGraw Hill Education New York ; 1995 . p. 154 – 161 . 46. ↵ Tu J , Yeoh GH , Liu C , Tao Y. Computational Fluid Dynamics: A Practical Approach . Elsevier ; 2023 . 47. ↵ Kuhn A , Fischer SC . On-lattice Vicsek model in confined geometries . arXiv . 2022 ; p. 2105.08792. doi: 10.48550/arXiv.2105.08792 . OpenUrl CrossRef 48. ↵ Obishakin E , Stijlemans B , Santi-Rocca J , Vandenberghe I , Devreese B , Muldermans S , et al. Generation of a Nanobody Targeting the Paraflagellar Rod Protein of Trypanosomes . PloS one . 2014 ; 9 . doi: 10.1371/journal.pone.0115893 . OpenUrl CrossRef 49. ↵ Bezanson J , Edelman A , Karpinski S , Shah VB . Julia: A fresh approach to numerical computing . SIAM Review . 2017 ; 59 ( 1 ): 65 – 98 . doi: 10.1137/141000671 . OpenUrl CrossRef PubMed 50. ↵ Danisch S , Krumbiegel J. Makie.jl: Flexible high-performance data visualization for Julia . Journal of Open Source Software . 2021 ; 6 ( 65 ): 3349 . doi: 10.21105/joss.03349 . OpenUrl CrossRef 51. ↵ Koirala RP , Dawanse S , Pantha N. Diffusion of glucose in water: A molecular dynamics study . Journal of Molecular Liquids . 2022 ; 345 : 117826 . doi: 10.1016/j.molliq.2021.117826 . OpenUrl CrossRef 52. ↵ Bowen WJ , Martin HL . The diffusion of adenosine triphosphate through aqueous solutions . Archives of Biochemistry and Biophysics . 1964 ; 107 ( 1 ): 30 – 36 . doi: 10.1016/0003-9861(64)90265-6 . OpenUrl CrossRef PubMed Web of Science 53. ↵ Chen C , Nakamura T , Koutalos Y. Cyclic AMP Diffusion Coefficient in Frog Olfactory Cilia . Biophysical Journal . 1999 ; 76 ( 5 ): 2861 – 2867 . doi: 10.1016/S0006-3495(99)77440-0 . OpenUrl CrossRef PubMed Web of Science 54. ↵ Zhang H , Freitas D , Kim HS , Fabijanic K , Li Z , Chen H , et al. Identification of distinct nanoparticles and subsets of extracellular vesicles by asymmetric flow field-flow fractionation . Nature Cell Biology . 2018 ; 20 ( 3 ): 332 – 343 . doi: 10.1038/s41556-018-0040-4 . OpenUrl CrossRef PubMed 55. Skliar M , Chernyshev V. Imaging of Extracellular Vesicles by Atomic Force Microscopy . Journal of Visualized Experiments . 2019 ; doi: 10.3791/59254 . OpenUrl CrossRef 56. ↵ Zhang P , Jiang J , Zhou X , Kolay J , Wang R , Wan Z , et al. Label-free imaging and biomarker analysis of exosomes with plasmonic scattering microscopy . Chemical Science . 2022 ; 13 ( 43 ): 12760 – 12768 . doi: 10.1039/D2SC05191E . OpenUrl CrossRef PubMed 57. ↵ Tronnolone H , Tam A , Szenczi Z , Green JEF , Balasuriya S , Tek EL , et al. Diffusion-Limited Growth of Microbial Colonies . Scientific Reports . 2018 ; 8 ( 1 ): 5992 . doi: 10.1038/s41598-018-23649-z . OpenUrl CrossRef PubMed 58. ↵ Slade AL , Cremers AE , Thomas HC . The Obstruction Effect in the Self-Diffusion Coefficients of Sodium and Cesium in Agar Gels . The Journal of Physical Chemistry . 1966 ; 70 ( 9 ): 2840 – 2844 . doi: 10.1021/j100881a020 . OpenUrl CrossRef 59. ↵ Ben-Jacob E , Cohen I , Gutnick DL . Cooperative organization of bacterial colonies: From Genotype to Morphotype . Annual Review of Microbiology . 1998 ; 52 (Volume 52, 1998): 779 – 806 . doi: 10.1146/annurev.micro.52.1.779 . OpenUrl CrossRef PubMed Web of Science 60. ↵ Armstrong JK , Wenby RB , Meiselman HJ , Fisher TC . The Hydrodynamic Radii of Macromolecules and Their Effect on Red Blood Cell Aggregation . Biophysical Journal . 2004 ; 87 ( 6 ): 4259 . doi: 10.1529/biophysj.104.047746 . OpenUrl CrossRef PubMed Web of Science 61. ↵ Moin P. Fundamentals of Engineering Numerical Analysis . Cambridge University Press ; 2010 . 62. ↵ Batchelor MT , Henry BI . Limits to Eden growth in two and three dimensions . Physics Letters A . 1991 ; 157 ( 4 ): 229 – 236 . doi: 10.1016/0375-9601(91)90057-F . OpenUrl CrossRef 63. ↵ Bogoyavlenskiy VA . How to grow isotropic on-lattice diffusion-limited aggregates . Journal of Physics A: Mathematical and General . 2002 ; 35 ( 11 ): 2533 – 2538 . doi: 10.1088/0305-4470/35/11/301 . OpenUrl CrossRef 64. ↵ Alves SG , Jr . SCF . Is it really possible to grow isotropic on-lattice diffusion-limited aggregates? Journal of Physics A: Mathematical and General . 2006 ; 39 ( 12 ): 2843 – 2852 . doi: 10.1088/0305-4470/39/12/001 . OpenUrl CrossRef 65. ↵ Paiva LR , Ferreira SC . Universality class of isotropic on-lattice Eden clusters . Journal of Physics A: Mathematical and Theoretical . 2006 ; 40 ( 1 ): F43 . doi: 10.1088/1751-8113/40/1/F05 . OpenUrl CrossRef 66. ↵ Yates CA , Baker RE . Isotropic model for cluster growth on a regular lattice . Physical Review E . 2013 ; 88 ( 2 ): 023304 . doi: 10.1103/PhysRevE.88.023304 . OpenUrl CrossRef 67. ↵ Verhoef B , Hermsen R , De Graaf J. Fluid-derived lattices for unbiased modeling of bacterial colony growth ; 2024 . View the discussion thread. Back to top Previous Next Posted November 05, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following An agent-based model of Trypanosoma brucei social motility to explore determinants of colony pattern formation Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share An agent-based model of Trypanosoma brucei social motility to explore determinants of colony pattern formation Andreas Kuhn , Timothy Krüger , Markus Engstler , Sabine C. Fischer bioRxiv 2025.11.04.686461; doi: https://doi.org/10.1101/2025.11.04.686461 Share This Article: Copy Citation Tools An agent-based model of Trypanosoma brucei social motility to explore determinants of colony pattern formation Andreas Kuhn , Timothy Krüger , Markus Engstler , Sabine C. Fischer bioRxiv 2025.11.04.686461; doi: https://doi.org/10.1101/2025.11.04.686461 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Biophysics Subject Areas All Articles Animal Behavior and Cognition (7637) Biochemistry (17705) Bioengineering (13899) Bioinformatics (41970) Biophysics (21463) Cancer Biology (18605) Cell Biology (25526) Clinical Trials (138) Developmental Biology (13385) Ecology (19911) Epidemiology (2067) Evolutionary Biology (24329) Genetics (15615) Genomics (22514) Immunology (17743) Microbiology (40424) Molecular Biology (17194) Neuroscience (88650) Paleontology (667) Pathology (2835) Pharmacology and Toxicology (4827) Physiology (7648) Plant Biology (15160) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)