Effects of Advective-Diffusive Transport of Multiple Chemoattractants on Motility of Engineered Chemosensory Particles in Fluidic Environments

Motility behavior of an engineered chemosensory particle (ECP) in fluidic environments is driven by its responses to chemical stimuli. One of the challenges to understanding such behaviors lies in tracking changes in chemical signal gradients of chemoattractants and ECP-fluid dynamics as the fluid is continuously disturbed by ECP motion. To address this challenge, we introduce a new multiscale numerical model to simulate chemotactic swimming of an ECP in confined fluidic environments by accounting for motility-induced disturbances in spatiotemporal chemoattractant distributions. The model accommodates advective-diffusive transport of unmixed chemoattractants, ECP-fluid hydrodynamics at the ECP-fluid interface, and spatiotemporal disturbances in the chemoattractant concentrations due to particle motion. Demonstrative simulations are presented with an ECP, mimicking Escherichia coli (E. coli) chemotaxis, released into initially quiescent fluids with different source configurations of the chemoattractants N-methyl-L-aspartate and L-serine. Simulations demonstrate that initial distributions and temporal evolution of chemoattractants and their release modes (instantaneous vs. continuous, point source vs. distributed) dictate time histories of chemotactic motility of an ECP. Chemotactic motility is shown to be largely determined by spatiotemporal variation in chemoattractant concentration gradients due to transient disturbances imposed by ECP-fluid hydrodynamics, an observation not captured in previous numerical studies that relied on static chemoattractant concentration fields.


Introduction
Intellectual and technological advances in a variety of fields continue to refine our understanding of the principles and potential applications of nanorobotic systems. Of great interest in this area is the understanding and development of control systems through which nanorobotic devices or bacterial biohybrids carrying a payload can be effectively directed to a specified target. Work in this field holds particular promise in applications relevant to health and biomedicine [1][2][3][4][5]. In addition to considerations of cargo and payload storage design, successful nanorobotics-based therapies must rely on a clear understanding of the structure and function of actuators, sensors, and power sources that govern the devices' abilities to interact with their fluidic environments [6]. While the fluidic environment poses navigational challenges to device design in each of these areas [7], natural biological systems, which have undergone adaptation and evolutionary selection for optimized solutions to these issues, have provided insights and inspiration [8][9][10][11][12]. Indeed, bacterial cells have been engineered to target specific locations in animal systems, most notably cancer tissue [2,[13][14][15].
A variety of actuation systems have been explored, including magnetic and acoustic fields, light and chemical energy [5]. While each of these actuating systems is characterized by distinct design features, advantages and limitations, they all share in common the fact that the engineered device operates in a fluidic environment. Therefore, it is essential to computationally and experimentally investigate how various fluidic environments influence the motility of particles, with the guiding principle that an understanding of actuation mechanics as well as particle-fluid dynamics will provide valuable insights into the design properties and behaviors of engineered devices in fluids. As such, biomimetics of the flagellar chemosensing bacterium Escherichia coli (E. coli), whose motility in fluids is driven by sensing chemical gradients and transformation of electrochemical energy into motion [16], could lead to the enhanced design of autonomous chemosensory structures with diverse applications in a variety of fluidic environments.
In this paper, we present a novel multiscale numerical model to simulate chemotactic behavior of an engineered chemosensory particle (ECP) swimming in a fluidic environment ( Figure 1). The model is built through the dynamic coupling of three modular components involving (i) ECP chemical attractant signal sensing and swimming response based on the chemotaxis adaptive phosphorelay circuit that governs the operation of the E. coli flagellar motor (MRC module); (ii) chemotactic ECP hydrodynamic interactions with the bulk fluid (CLB module); and (iii) chemical substrate transport phenomena in a fluidic environment, validated against a 2D benchmark problem and described in this paper (ADT module). The first module (MRC) simulates submicron−scale cell signaling processes that determine the run and tumble biased random walk behavior of the ECP based on the chemical environment. The second and third modules (CLB and ADT) simulate particle-fluid hydrodynamics and spatiotemporal variations in the fluid velocity and substrate concentration at the cm-scale. Hence, the new model is a coupled multiscale numerical model that simulates the transformation of chemical energy by an ECP to mechanical energy as it swims through fluidic environments containing concentration gradients of chemoattractants. This is accomplished by accommodating dynamic changes in spatiotemporal distributions in the fluid velocity and concentration fields due to ECP motion. Details of each of these module components comprising the multiscale model are described below.
We recently developed multiscale computational models for E. coli chemotactic sensing by coupling the Monod-Wyman-Changeux mixed chemosensory receptor cluster model, known as RapidCell (RC) [17], with the method of regularized Stokeslets [18] or the colloidal lattice Boltzmann (CLB) models [19,20]. These models were used to simulate motility of chemosensory particles in confined fluidic environments with externally imposed chemoattractant gradients [21,22]. In this paper, we extend the coupled RC-CLB model to simultaneously simulate advective-diffusive transport (ADT) of two unmixed chemoattractants, ECP-fluid hydrodynamics, and disturbances in spatiotemporal distributions of chemoattractants due to particle motion. We equip the ECP with the well-established model of chemical stimulus sensing circuitry of E. coli, which confers on the ECP the ability to detect and respond to gradients of chemoattractant compounds in the fluidic environment. We refer to our upgraded model as MRC-CLB-ADT, corresponding to the modified RapidCell (MRC)-colloidal lattice Boltzmann (CLB)-advective-diffusive chemoattractant transport (ADT) model. The MRC-CLB-ADT model is used in this paper to follow the position and swimming trajectories of an ECP in fluidic environments in which two unmixed chemoattractants are introduced on opposite sides of a confined flow domain. We also utilize the model to investigate how the residence times of the ECP are affected by chemoattractant release mode (instantaneous vs. continuous, point source vs. distributed) as concentration and fluid velocity fields are altered through ECP swimming behavior.
In the following sections, we first describe the mathematical formulation of the new coupled MRC-CLB-ADT model to simulate the motility of ECPs in two-dimensional (2D) fluidic environments in response to spatiotemporal variations in concentrations of the amino acid chemoattractants N-methyl-L-aspartate (MeAsp) and L-serine (Ser). The selection of these amino acids is based on their frequent use in experimental studies of E. coli chemotaxis mechanics and chemoreceptor biochemistry [16]. Using this modeling framework, we present and discuss the results from simulations with increasing levels of complexity and realism. The first set of simulations is performed using imposed temporally-invariant, but spatially-variant concentration fields. In these simulations, chemoattractant distributions are presumably not affected by ECP-fluid hydrodynamics, and are used as benchmark cases. A second more realistic set of simulations incorporates the effects of ECP motility in the fluidic environment on the chemoattractant distributions. Through the incorporation of a dynamically-evolving mixed chemical environment, the MRC-CLB-ADT model simulations reveal critical roles of the ECP-fluid hydrodynamics on the chemosensory particle motility not previously recognized in static concentration field-based models.

Mathematical Framework of the MRC-CLB-ADT Model
A mathematical description of main submodels (MRC, CLB, and ADT) and the optional submodel (static chemoattractant concentration fields) of the MRC-CLB-ADT model and their coupling are provided in this section. Numerical validations of the CLB and ADT modules are also discussed in this section.

Module 1. Modified RapidCell (MRC) Model for Particle Chemosensing in Two Chemoattractant Fields
We modified the RapidCell (RC) model to simulate chemotactic motility of the ECP with the premise that the ECP mimics E. coli chemotaxis in the presence of two unmixed chemoattractants. The RC model [17] was originally developed to simulate flagellar bacterial chemotaxis in an environment with a spatiotemporally varying concentration gradient of a single chemoattractant, and performs two tasks: chemoattractant signal processing by the methyl-accepting chemotaxis protein (MCP) sensory lattice and adaptive feedback response of the sensory array [16,[23][24][25]. Signal processing by the cell occurs through interactions between chemoattractant-activated chemoreceptors, CheA kinase, and other regulators including CheY, CheR, CheB, and CheZ ( Figure 2). Chemoattractant-receptor interactions regulate the autophosphorylation activity of CheA [24]. CheA-P phosphoryl transfer activity may then result in switching of direction of motor rotation through CheY-P [26] or adjustment of adaptive response through MCP methyl group hydrolysis by CheB-P [27]. coli. Chemoattractants such as N-methyl-L-aspartate and L-serine (represented by triangle and diamond shapes) are sensed by Tar and Tsr MCPs, respectively, and binding results in signal transduction across the cell membrane to a phosphorelay response circuit. Phosphoryl group (P) transfer to Che proteins controls direction of rotation of flagellar motor and MCP methylation-dependent adaptive response. Default flagellar rotation is counterclockwise, causing cell to run; switching to clockwise rotation results in reorientation of cell through tumbling motion due to flagellar unbundling. Circled letters represent Che proteins described in text, e.g., A represents CheA, Y represents CheY, etc.
The RC model was modified to account for the presence of two unmixed chemoattractants using the total free energy differences in [28]. The effect of the total free energy differences between 'on' and 'off' states for two receptors sensing two chemoeffectors is described as where N is the number of chemoreceptors in the receptor cluster.
[MeAsp] and [Ser] are the extracellular chemoattractant MeAsp and Ser concentrations, respectively [28]. The binding of MeAsp by Tar is given in the modified total free energy by ν a and the binding of Ser by Tsr is represented in the modified total free energy by ν s . The offset energy, h(m) is given by 1 − m/2 where m is receptor methylation defined in Equation (4). The dissociation constant for the chemoattractant in the 'on' or 'off' state is specified as K on/o f f r (r = a, s). The total free energy differences F from Equation (1) are used in the RC model to compute the receptor methylation (m) and the basal motor bias (mb) where A c is the probability of the cluster activity, CheY-P is the concentration of phosphorylated CheY, k Y , k Z and γ Y are the rate constants, k s is the scaling coefficient, [CheZ] tot is the total CheZ concentration, a and b are the methylation scaling factors, m is receptor methylation, mb is the motor bias, mb 0 is the basal motor bias, and H is the motor Hill coefficient. In RC model simulations, the initial methylation state of the receptor cluster is obtained from a steady-state methylation level associated with the initial chemoattractant concentration bound to the chemoreceptor cluster. The methylation is then updated by solving Equation (4) using the forward Euler finite difference method. Motor bias values may range from 0 to 1 with higher values corresponding to a greater likelihood of running motion of bacteria. To determine whether a particle will run or tumble, a uniform random number ξ is generated between 0 and 1; and if ξ < mb, the particle runs; otherwise, it tumbles. A complete list and description of the parameters, variables, and functions in Module 2.1 are provided in Tables A1 and A2 in the Appendix A.4.1.
Our modification to the RC model provides a framework for relating the time history of multiple chemoattractant concentration sensing events at the particle's chemosensory array to the run-tumble probability output and is called the modified RapidCell model (MRC).

Static (Time-Invariant) Concentration Fields
If externally-computed time-invariant radial concentration gradients are used for MeAsp and Ser in numerical simulations, the concentration gradients are described by [28] [MeAsp] = ω C a0 where C a0 and C s0 are the minimum chemoattractant concentrations for MeAsp and Ser, respectively.
x and y are the horizontal and vertical coordinates, and ω and ν are scaling parameters for MeAsp and Ser gradients, respectively. The parameters (x a , y a ) and (x s , y s ) describe the location of the maximum MeAsp and Ser concentrations in a 2D square domain, described as (x, y) = {x, y ∈ R : [−L * , L * ]}, in which L * is the domain length [28]. The imposed chemoattractant concentration gradients of MeAsp and Ser are scaled, using a scaling parameter of r, for larger domains in RC-CLB simulations as follows [Ser] = νr C s0 + exp − (x + x s ) 2 + (y + y s ) 2 r .

Dynamic (Time-Variant) Concentration Fields
In reality, the MeAsp and Ser concentrations change according to their own diffusion process while being advected by the fluid flow and the ECP motion as time goes on. Therefore, the concentration fields cannot be externally computed and independent of time as in the ideal case of Section 2.1.1. At each time step the fluid-ECP interactions need to be taken into account so that the distributions of MeAsp and Ser concentrations can be updated properly. Section 2.2 presents how the fluid-ECP interactions are modeled. Section 2.3 shows how the resultant fluid velocities from Section 2.2 affect the transport of MeAsp and Ser.

Module 2. Colloidal Lattice Boltzmann (CLB) Model for Particle-Fluid Interactions
The CLB model has two submodules, (i) fluid flow submodule and (ii) particle flow submodule, The former calculates the local changes in the fluid velocity field in response to ECP motion. The latter subsequently updates angular and translational velocities of an ECP in the disturbed fluid velocity field. These two new submodules operate sequentially in each time step and calculations are based on momentum exchanges between a motile ECP and the bulk fluid across the ECP's surface. In the fluid flow module, the fluid flow is governed by the Navier-Stokes equation and it is solved using the lattice-Boltzmann method (LBM) [29]. In the particle flow module, the angular and translational velocities of an ECP are computed based on Newton's equation of motion. Mathematical details of these two submodels are presented next.

Fluid Flow Submodule (FFS)
In the fluid flow submodule of the CLB model, the mesodynamics of the transient, weakly compressible, Newtonian fluid flow is described by a single relaxation time [30].
where f i (r, t) is the complete set of population density of discrete velocities e i at position r and discrete time t with a time increment of t, and τ is the relaxation parameter. The left-hand side of Equation (8) describes the streaming of populations from a lattice node r to the closest neighbor in the direction e i . The right-hand side of Equation (8) represents the local collision process. f eq i in Equation (8) is the discrete equilibrium Maxwell-Boltzmann distribution approximated by the low Mach number mass and momentum conserving expansion [31] where ω i is the weight associated with e i , c s is the speed of sound, c s = x/ √ 3 t, and x is the lattice spacing. The local fluid density, ρ, and velocity, u, at the lattice nodes are given by where g is the acceleration due to external forces [32]. In Equation (10), ρg = P/L * , and hence, ρg specifies the pressure differential across the fludic domain with the length of L * . If the fluidic domain is stagnant, then g = 0. Otherwise, the flow strength across the fluidic domain can be specified by ρg.
Through the Chapman-Enskog approach [29], in the limit of small Knudsen number for weakly compressible fluids ( ρ/ρ ∼ M 2 ∼ 1 × 10 −4 , where M is the Mach number), the LB method for single-phase fluid flow recovers the Navier-Stokes equations whereν is the kinematic viscosity of the fluid. Pressure is computed via the equation of state for an ideal gas, P = c 2 s ρ. A commonly used D2Q9 (two-dimensional nine velocity vector) lattice geometry (Figure 3), which ensures fourth-order lattice isometry, was adopted in LB simulations in this paper. For which unit discrete velocities e i at each lattice node are defined as The first column vector of e in Equation (12) is the null vector corresponding to the rest population, the second through fifth column vectors correspond to the four vectors of length unity directed toward the nearest neighboring lattice nodes, and the sixth through ninth column vectors correspond to the four vectors of length √ 2 directed toward the next-nearest neighboring lattice nodes ( Figure 3). Equation (8) recovers the Navier-Stokes equation for ω i = 4/9 for the rest populations (i = 0), ω i = 1/9 for the off-diagonal populations (i ∈ [1,4]), and ω i = 1/36 and for the diagonal populations (i ∈ [5,8]) in Equation (9) and τ = 0.5 + 3ν t/( x) 2 in Equation (8), which are obtained through the Chapman-Enskog expansion [29]. Hence, τ in the LB method is determined by the kinematic viscosity of the fluid,ν. Numerical stability in fluid flow simulations was ensured by τ = 0.8 < 1.0.
In each time-step in numerical simulations, f i and f eq i were computed at each lattice node via Equations (8) and (9). f i 's can be altered locally by ECP motion, which will be discussed in the next section. After f i 's were computed, the fluid velocity and density at each lattice node were computed via Equation (10). In these simulations, the fluid was Newtonian; therefore,ν, and hence, τ were constants. Moreover, x and t remained constant throughout simulations.
The LB method is second accurate in space and time. The LB model was preferred over classic Navier-Stokes solvers in this paper as the computationally demanding nonlinear convective term in the Navier-Stokes equation, (u · ∇) u, is replaced by a linear arithmetic streaming term (the left-hand side of Equation (8)) in the LB method. The streaming is exact and local mass and momentum conservations are accurate to machine round-off error [29,33].

Particle Flow Submodule
The particle flow submodule of the CLB model [19,20], built on the LBM formulation by [34][35][36], was modified to simulate hydrodynamic interactions between an ECP (represented as a circular-cylindrical particle in the LBM) and the bulk fluid [21]. ECP-fluid hydrodynamic forces, F r b at boundary nodes located halfway between the intra-particle lattice node, r v , and extra-particle lattice node, r v + e i , are computed based on momentum exchanges between the ECP and surrounding fluid ( Figure 4) [34,37] where f i is the population density in the −e i direction at the post-collision time t * , and u r b is the local particle velocity at the boundary node r b . Equation (13) is related to the continuum-scale particle-fluid hydrodynamic force on the particle where m p is the particle (ECP) mass and U p is the particle velocity. Steric interaction forces, F r i , between the particle and stationary solid zones are expressed in terms of two-body Lennard-Jones potentials [38] such that n, where | r i | is the distance between a particle surface node and the stationary solid node located on channel walls or inline obstacles (r i = r pw ); p is the ECP index; w is the wall or obstacles index; | r it | is the repulsive threshold distance; n is the unit vector along r i ; and ψ is the stiffness parameter used to adjust the repulsive strength between the particle and stationary solid zones. Then, the total force, F T , on the ECP is , and uncovered, r u b , lattice nodes due to particle motion [35,37,39]. The force associated with the running motion of the ECP is defined as F run = (r cl −r c ) |r cl −r c | f m , where r cl is the location of the receptor cluster, r c is the position of the particle centroid, and f m is the force strength. r cl is placed on the boundary nodes. The total torque on the ECP, T T , is defined as  [34,38,39]. Filled circles are the intra-particle virtual fluid nodes of ECP closest to its surface, filled triangles outside the ECP surface represent its extra-particle bulk fluid nodes, and the filled square represents the boundary node located half-way between the intra-particle node (r v ) and extra-particle node ( r v + e i t). Hydrodynamic links along which the ECP exchanges momentum with the fluid are shown by red line segments.
Torque due to the ECP tumbling motion as it seeks the highest chemoattractant gradients, T tumble , is where Υ is the time-scaling factor associated with the angular rotation of the ECP due to its tumbling motion. This is different from the time-scale associated with its translation velocity resulting from running motion. Ω tumble (t) is the ECP's angular velocity due to particle torque resulting from its tumbling state at time t and is defined as Ω tumble (t + ∆t) = Ω tumble (t) + tT tumble (t + t) /I p , where I p is the moment of inertia of the ECP (I p = m p r p 2 /2, in which r p is the radius of the circular ECP) and Ω tumble (t = 0) = 0. The angular rotation, θ, of the ECP associated with its tumbling motion over t is computed by θ = 2π(ϕ − 0.5) where ϕ is a uniformly distributed random number between 0 and 1. The translational velocity, U p , and the angular velocity of the ECP, Ω p , are advanced in time according to the discretized forms of Newton's equations of motion Local velocities at the boundary nodes are computed by The new position of the ECP is computed by Then the position of the receptor cluster on the ECP surface is updated by where r p is the particle radius and the rotational angle of the receptor cluster is defined as Finally, the population densities at the intra-particle node, r v , and the extra-particle node, r v + e i t, are updated to account for momentum-exchange between the ECP and bulk fluid via [34] ( Figure 4) In each time-step in numerical simulations, the total force, F T , and the total torque, T T , on the ECP were computed via Equations (15) and (16). The resultant translational and angular velocities of the ECP were calculated by Equations (18) and (19). Local velocities at the boundary nodes of the ECP and the new position of the ECP were computed next by Equations (20) and (21). The new location of the receptor cluster on the ECP's surface, at which chemical sensing of chemoattractants was simulated via MRC, was then computed by Equation (22). Local disturbances in the immediate vicinity of the motile ECP altered population densities, f i , in accordance with Equations (23) and (24). The altered f i 's near the ECP's surface were used to compute the new fluid velocity field via Equations (8) and (10).
The CLB model was validated in [36] against numerically-computed (via the finite-element method) settling trajectories of a circular-cylindrical particle in an initially quiescent fluid [40]. As reported in [36], the CLB model also closely predicted terminal velocities of spherical particles 5% and 10% denser than the bulk fluid [36] in particle settling experiments reported in [41].

Module 3. Advective-Diffusive Transport (ADT) Model for Chemoattractant Distributions
In the original RC model [17], the chemoattractant environment surrounding chemosensory particles was not fluid-based; therefore, static chemottractant concentrations were externally computed (as discussed in Section 2.1.1) and artificially imposed onto the fluidic domain. To overcome this shortcoming, a chemoattractant transport model, based on the LBM, was formulated in this paper to simulate spatiotemporal distributions of chemoattractants in the fluidic environment of the ECP by accommodating the effects of particle motion-induced disturbances in the flow and chemoattractant concentration fields. The advective-diffusive chemoattractant transport was solved using the LBM on a D2Q9 lattice [42,43]: where g i (r, t) is the complete set of population density of discrete velocities e i associated with the chemoattractant concentration at position r and time t with a time increment of t. τ c is the relaxation parameter associated with the chemoattractant transport. In Equation (25), g eq i is the local equilibrium for the chemoattractant transport process and is given by [43] The local chemoattractant concentration at each lattice node r is given by C = ∑ i g i . Through the Chapman-Enskog expansion [43], Equations (25) and (26), recovers the continuum-scale transient advective-diffusive substrate transport equation, where C is the chemoattractant concentration, D is the Fickian diffusion coefficient, and u is computed by Equation (10). Equation (25) is equivalent to Equation (27) The ADT model was validated with a 2D benchmark problem (Appendix A.1) and its performance was tested with different flow simulations (Appendix A.2) in the Appendix A.
In each time-step in numerical simulations, g i and g eq i were computed at each lattice node via Equations (25) and (26), using the fluid velocity computed at each lattice node by Equation (10). After g i 's were computed, chemoattractant concentrations at each lattice node were calculated by C = ∑ i g i . Computed chemoattractant concentrations were used by the MRC to determine the tumble or run motion of the ECP. Decision on the tumble or run motion of the ECP affected F T and T T in Equations (15) and (16), and hence, altered local f i 's in the vicinity of the mobile ECP according to Equations (23) and (24). A complete list and description of variables and parameters used in the CLB and ADT modules are provided in Table A3.
In this paper, we adopted the Fickian diffusion process via Equation (27), following the earlier numerical work by [42,43], in which D is constant in space and time, so that the diffusion process is described by D∇ 2 C. However, future work, involving also experimental tasks, will explore the use of the Fokker-Planck equation instead to describe the diffusive processes with temporal dispersion rather than spatial dispersion [44][45][46][47] of chemoattractant transport in determining chemotactic motility of ECPs. The LB method has been shown to be capable of solving 2D Fokker-Planck equations with variable coefficients by [48].

Coupling of the Modules, MRC-CLB-ADT Model
The MRC and CLB models are coupled via Equations (15), (16), and (22) [21]. Although each ECP is subject to its own signal processing and chemotactic swimming behavior, as described by Equations (1)-(5), their interactions with the bulk fluid and stationary solid zones are accounted for by particle-fluid hydrodynamic forces and particle-wall steric interaction forces in Equation (15). If the concentration fields for MeAsp and Ser are externally computed and imposed onto the fluidic domain using Equations (6) and (7), the coupled MRC-CLB model may be used to describe the transient behavior of ECPs without using the ADT model.

Simulation Parameters
The motility of chemosensory particles in a fluidic environment is typically described by Reynolds number, Re = |U p |D p /ν where D p is the particle size. In our simulations, considering the average |U p | being 6.09 × 10 −3 cm/s,ν being 8.70 × 10 −3 cm 2 /s, and the representative size of the ECP being 2 cm, the Re associated with the ECP motility is 1.4, residing in the creeping flow regime. In this case, the ECP travels from one side of the square flow domain of 3232.4 cm 2 to the opposite side in a straight flow path in 2.6 h at a rate of 0.3% of its body length per second.
ECP's mass was 4.2 g in simulations. The force strength associated with the run motion of ECP was set to 3.25 × 10 −4 kN to keep Re on the order of 1. The time-scaling factor associated with tumbling-induced rotation of bacteria, Υ, in Equation (17), was previously estimated to be 5 from simulations with a single chemotactic particle released into a confined domain in the absence of inline obstacles [21]. Therefore, Υ = 5 was adopted in simulations in this paper.
In simulations, |r pw | = 1.5 lattice unit (l.u.) and ψ = 1. ECP-wall steric interaction forces would be non-zero only when surface boundary nodes of the ECP move within 1.5 l.u. (0.43 cm) of stationary wall surfaces to avoid physically unrealistic overlaps [36]. When the boundary nodes of ECP the moves in within 1.5 l.u. of wall nodes, instantaneous, short-lived (within a few time-steps) relatively large (typically within an order of magnitude of F r b ) steric pulse applies to keep the ECP-wall separation distance larger than 1.5 l.u.
Our simulations were conducted with two unmixed chemoattractants, MeAsp and Ser. When these concentrations were imposed in the environment using the Equations (6) and (7), the scaling parameters ω and ν controlled the peak concentrations. ω = 1 returned the optimal sensitivity of chemoreceptors to MeAsp in the imposed environment as observed in [28]. Two values of the scaling parameter ν (0.1 or 0.001) were chosen to control the sensitivity of the chemoreceptor response to Ser. When ν = 0.1, a dense level of cell accumulation to the Ser peak was expected in [28] due to a higher peak concentration compared to the one with ν = 0.001. Therefore, in our simulations, ω was kept at 1 while ν = 0.1 or ν = 0.001 to compare the transient behavior of ECPs in different scenarios given the imposed environment (i.e., Equations (6) and (7) were used instead of the ADT model) or the dynamic environment (i.e., the ADT module was coupled as shown in Figure 5).
In our simulations, the diffusion coefficient values, D, for Ser and MeAsp were chosen such that they are within the range of experimentally reported values for Tar and Tsr substrates [49][50][51]. Because diffusion of chemoattractants would be enhanced by disturbances (additional mixing) in the fluid [52] by ECP motion, the values of D were increased by a factor of 10 to account for enhanced diffusion in MRC-CLB-ADT simulations. The factor of 10 was chosen so that the concentration distributions were spreading out reasonably within 50,000 time steps. However, the actual value of the enhancement factor and D could be obtained from experiments in future studies. Thus, for our simulations, D for Ser is set to 8.7 × 10 −5 cm 2 /s and D for MeAsp is set to be 1.08 times higher than for Ser.
The scaling parameter r in Equations (6) and (7), is set to 14.3 cm. The initial maximum concentrations of MeAsp and Ser are located at (14.3 cm, 28.9 cm) and (42.9 cm, 28.9 cm), respectively, which correspond to P 1 = (50, 101) and P 2 = (150, 101) on a lattice grid. The minimum concentrations for MeAsp and Ser, C a0 and C s0 , are set to 0.1 µM. In our simulations, 1 unit of concentration corresponds to 1 µM.

Results
In MRC-CLB-ADT simulations discussed in the subsequent sections, except for the validation test

Simulations with Imposed Temporally-Invariant, Spatially-Variant Chemoattractant Concentrations
Base case simulations involving temporally-invariant, spatially-variant chemoattractant concentration fields were used to compare MRC-CLB model results to the simulation results by Edginton and Tindall [28]. Time-invariant chemoattractant concentration fields ( Figure 6) computed by Equations (6) and (7) were imposed onto the flow field, instead of being computed in each time-step by the ADT model. Hence, temporal disturbances in the chemoattractant concentration fields due to motility of the ECP were not accounted for. These base case simulations are classified as "Imposed" due to the static nature of the chemoattractant concentration fields artificially imposed onto the fluid domain.
The MeAsp gradient parameter was set to ω = 1, based on previously reported Tar sensitivity curve data in [54]. The gradient parameter for Ser was set to ν = 0.1 or ν = 0.001. These ν values, in combination with the fixed ω gradient value, were chosen to account for changes in the bias of an ECP to travel toward increasing MeAsp concentrations due to saturation effects [28,54]. To account for the biased motion of an ECP in this environment, ten trial simulations were performed with ν = 0.1 or ν = 0.001. The ECP trajectories were different in these replicates due to the randomness in the MRC method (Equation (17) (101, 101).
The coordinates (x, y) of the centroid of a motile ECP were tracked for 50,000 time steps, corresponding to 13 h. At any given time-step, if x > 100, the ECP would be located in the right half of the domain, where the center of the Ser concentration field was initialized at (150, 101), which will be referred to as "On Ser Half". Similarly, if x < 100, the ECP would be located in the left half of the domain, where the center of the MeAsp concentration field was initialized at (50, 101), which will be referred to as "On MeAsp Half" hereafter.  Figure 7 shows the average number of time steps from ten replicate simulations, for which the ECP was found in either the Ser-or MeAsp-half of the domain, in response to static chemoattractant gradients specified by ν = 0.1 or ν = 0.001. The large error bars indicate that trajectories of the ECP from these replicates were broadly distributed across the domain. The overall trends of the average values and error bars in Figure 7 are in good agreement with the results by [28], which showed an affinity of the chemosensory particle to (i) the Ser concentration field with scaling parameters ω = 1 and ν = 0.1, and (ii) the MeAsp concentration field with scaling parameters ω = 1 and ν = 0.001. These results reveal that ECP motility is governed by receptor sensitivity rather than absolute chemoattractant concentrations when time-invariant concentration fields are assumed.

Simulations with Spatiotemporal Variations in Chemoattractant Concentrations Computed via ADT Model
Using the coupled MRC-CLB-ADT model ( Figure 5), two-dimensional simulations of ECP motility with incrementally improved realism in the problem set-up are discussed in this section. In order to represent the chemotatic behavior of an ECP in fluidic environments more realistically, the ADT model is used to simulate spatiotemporal variations in Ser and MeAsp concentration fields as the fluid is continuously disturbed by ECP motion. This is a significant conceptual and modeling improvement over the simulations discussed in Section 3.1, in which externally-computed time-invariant chemoattractant concentrations were artificially imposed onto the fluidic environment. Here, four cases with different chemoattractant source specifications were implemented in MRC-CLB-ADT simulations to analyze the effect of the initial distributions and release modes of chemoattractants on the transient chemotactic motility of an ECP:  • Case 2. "ADT Point: Continuous": After the initial condition was set up as in Case 1, ( C) of each chemoattractant was released into the fluid in each time-step for t > t from their respective point source locations, at which their maximum concentrations were maintained throughout the simulation. Snapshots from this simulation are shown in Figure 9. • Case 3. "ADT Imposed: Initial": Initial concentration fields of the chemoattractants were specified via Equations (6) and (7) and imposed onto the fluidic domain. No additional chemoattractant releases occurred for t > 0. Snapshots from this simulation are shown in Figure 10. • Case 4. "ADT Imposed: Continuous": After the initial concentration fields of chemoattractants were established as in Case 3, ( C) of each chemoattractant was released into the fluid in each time-step for t > t from their respective point source locations. Snapshots from this simulation are shown in Figure 11. In early times, the ECP moved toward the Ser field on the right half of the fluidic domain. However, as the Ser concentration gradually diffused out, the trajectory of the ECP became unpredictable. The motile ECP continuously disturbed and altered the flow field and concentration fields of the chemoattractants as it exchanged momentum with the bulk fluid. The resultant alterations in the concentration fields subsequently affected the signaling pathway of the ECP, and hence, its tumbling and running motion. These simulations demonstrate that trajectories of the ECP were sensitive not only to initial distributions of the chemoattractants, but also to temporal variations in the chemoattractant concentration fields. Therefore, the consideration of ECP motion-induced disturbances in the fluid velocity field is imperative in ECP motility studies, and hence, should be accommodated in simulations. Ensemble-average of temporal changes in chemotactic activities of the ECP over ten replicates for Case 1 through Case 4 are shown in Figure 12. The ECP was released from the center of the bounded domain and its position was tracked for 50,000 time steps. In Figure 12, the distance between the ECP's location at any given time-step and the center (50, 101) of the MeAsp field at t = 0 was computed and then averaged over ten replicates. The ensemble-average distance is denoted by d (50,101) . Similarly, the distance between the location of the ECP at any given time-step and the center (150, 101) of the Ser field at t = 0 was computed and then averaged over ten replicates. The ensemble-average distance is denoted by d (150,101) . Initially, the release locations or spatial distributions of the chemoattractants and the release locations of the ECPs were symmetric about the midpoint of the fluidic domain. Therefore, in each time-step, if [d (50,101) −d (150,101) ] < 0, the ECP would be on the left half of the domain, where the center of the MeAsp concentration field was initialized ("On MeAsp half"); otherwise it would be on the right half of the domain, where the center of the Ser concentration field was initialized ("On Ser half").
The "Imposed" case in Figure 12a exhibits consistent results with the the bar graph in Figure 7 in the sense that the ECP motility was biased toward the MeAsp/left half for ν = 0.001 (solid red curve), but toward the Ser/right half for ν = 0.1 (dashed blue curve). Figure 12a further reveals that in the end of the simulations, the ECP tends to remain on the side (either Ser or MeAsp) from which it was initially released. The early steep rise for ν = 0.1 is associated with the strong initial response of the Tsr chemosensory component to the static Ser gradient. The initial average ECP behavior in the "Imposed" condition with ν = 0.001 also shows a strong tendency toward the MeAsp gradient, albeit with a lower sensitivity toward this attractant. This could be due to the lower sensitivity of the Tar receptor relative to that of the Ser receptor Tsr. Intermittent periodicity in ECP's ensemble-average behavior in Figure 12a is associated with the sequence of its tumble and run motion. As the ECP continually runs and tumbles, it would move into, through, and out of zones of maximal chemoattractant concentration in some time steps. When the ECP senses continual temporal decreases in chemoattractant concentrations where it resides in the fluidic environment, it would begin to tumble and run back up to the zones with steeper concentration gradients. Based on the average positions of the ECP in reference to the center of the chemoattractant field in Figure 12a, such correction occurs sooner with Ser chemosensing for ν = 0.1 than MeAsp sensing with parameter ν = 0.001. Simulations with ν = 0.1 and ν = 0.001 showed in Figure 12b-d that the ECP in general traveled across both halves of the fluidic domain, rather than residing mostly in one half the domain as in Figure 12a, when spatiotemporal disturbances in the concentration fields due to ECP motion are accounted for. Simulations with Case 1 (Figure 12b) or Case 4 (Figure 12e) show that the ECP maintained its transient chemotactic activities predominantly on the MeAsp half, as in the "Imposed" Case, when ν = 0.001. On the other hand, the effect of this enforced biased travel on the ECP's ensemble-average trajectories was negligible in Case 2 (Figure 12c) The ECP maintained more than 99% of its chemotactic activities on the Ser half. Furthermore, the point injection of chemoattractants to the fluidic domain, in comparison to the imposed concentration fields at t = 0, resulted in transient chemotactic activities of the ECP on the opposite halves of the fluidic domain for ν = 0.1. Even in Case 3, where the initial concentration fields were set up as in the "Imposed" case, the concentration fields were spatially and temporally evolving due to ECP motion, which affected run and tumble motion and overall ensemble-average trajectories of ECP. Similarly, bias in ECP's motility toward the Ser half introduced by ν = 0.1 had insignificant effect on particle ensemble-average trajectories in Case 1 and Case 4 when the concentration fields are not static. As a result, the ECP retained most of its chemotactic activities in the MeAsp half. In contrast, the biased ECP traveled toward the Ser half as in the "Imposed" case.
In summary, Figure 12 demonstrates that the biased ECP motility in a fluidic environment with multiple static chemoattractant concentration fields (i.e., "Imposed" environment) can be enforced and controlled through a fixed decisive factor of ν. However, if the concentration fields dynamically evolve spatially and temporally, ECP motility and trajectories are determined by its interactions with dynamically-evolving surrounding fluidic and chemoattractant environments, which cannot be represented accurately by a static factor ν.

Discussion
The average number of time steps the ECP spent in the MeAsp half or Serine half (i.e., the mean residence time of the ECP in each half) computed by the MRC-CLB-ADT model for Case 1 through Case 4 is compared against the "Imposed" case in Figure 13. Although ECP trajectories in Case 3 and the "Imposed" simulations were different in Figure 12a,d, the resultant mean residence times of the ECP and the associated errors were comparable in Figure 13. As noted previously, Case 3 and the "Imposed" had the same initial chemoattractant distributions, but the concentration field dynamically evolved in Case 3, unlike in the "Imposed" simulation . Although trajectories of the ECP were different in Case 3 and the "Imposed' in each realization, the ensemble-average (the statistical mechanics) ECP motility data were similar and appear to be insensitive to spatiotemporal variations in the chemoattractant fields when the chemoattractants are spatially distributed initially and no additional chemoattractants are introduced into the fluid in later times. Hence, in such problem set-ups, the mathematically simpler and computationally less-expensive "Imposed" case could be used to evaluate the statistical mechanics of ECP motility in the design or performance assessment of ECPs.
When the chemoattractants were continuously released into the flow field in Case 4, the ECP preferentially spent more time in the MeAsp half, regardless of the value chosen for ν ( Figure 13). This rather surprising result can be explained by the concentration distributions in Figure 11. Case 4 simulation is similar to the "Imposed" simulation with the exception that C of the chemoattractants was continuously injected at the point sources in each time step. Although the injection maintained the highest local concentration of chemoattractants at the source locations, MeAsp spread more rapidly than Ser throughout the domain due to its higher diffusion coefficient, which in turn more strongly influenced the chemotactic activities of the ECP. This is clear from Figure 12e that shows a strong bias of ECP trajectory toward the MeAsp half over the entire simulation period.
When the chemoattractants were only initially released into the fluid as point sources in Case 1, the ECP spent more time in the MeAsp half of the domain. The maximum chemoattractant concentrations occurred at the point source location only at the start (Figure 10b). Because the ECP was initially farther away from point source locations and no additional chemoattractants were released into the fluid at later times, the ECP was incapable of accurately detecting the chemoattractants in early times. However, as advection and diffusion processes redistributed the chemoattractants as the time advanced, ECP trajectories were governed by continuously spreading and gradually diminishing local concentration gradients of chemoattractants.
(a) (b) Figure 13. Ten simulations of 50,000 time steps each were performed for the"Imposed" case, and Case 1 through Case 4, in which trajectories of an ECP in a fluidic environment with an ω/ν ratio of (a) 1/0.1 and (b) 1/0.001 were traced. Heights of the bars correspond to the total residence time of an ECP either in the right half, in which the Ser concentration was initialized, or in the left half, in which the MeAsp concentration was initialized, of the fluidic domain. The first set of bars in both (a) and (b) are repeated from Figure 7 while the remaining sets are from ADT model simulations.
When the chemoattractants were continuously released into the fluid for t > 0 in Case 2, the ECP retained its chemotactic activities mostly in the Ser half, regardless of the value chosen for ν. A comparison of Case 1 and Case 2 simulations with ν = 0.001, (Figure 10b,c) shows that the ECP in Case 2 was initially positioned too far from MeAsp to detect any signal and respond to it. As a result, despite the smaller Ser gradient parameter value, the ECP responded to Ser only. Apparently, Tsr sensitivity dominated over Tar sensing even for the lower ν value. Similarly, ECP response to Ser in simulations with ν = 0.001 was delayed as the ECP was initially too far from the Ser point source location (Figure10c, ν = 0.1). However, stronger Tsr response was observed after local chemoattractant gradients were established in the fluid and sensed by the ECP at approximately time step 15,000. These findings provide compelling evidence that unlike in Case 3, the "Imposed" condition would not be suitable to simulate chemotactic behaviors of the ECP in the problem set-ups in Case 1, Case 2, and Case 4.
In summary, unlike previously reported simulation results with imposed static chemoattractant concentration fields [28], our results with dynamically changing chemoattractant fields in response to ECP motion reveal for the first time that the relative magnitudes of ω and ν are not the sole factors in determining on which side of the fluidic domain an ECP would reside in. Instead, mutual dynamic interactions between particle motion and dynamically-varying concentration and flow fields would determine the statistical mechanics (ensemble-average) of ECP motility. Moreover, the release modes of the chemoattractants (point vs. non-point source and/or initial vs. continuous releases) and spatiotemporal evolution of chemoattractant concentration fields are found to be the critical factors for chemotactic activities and the statistical mechanics of motility of ECP, which would determine the zone(s) where an ECP would reside in a fludic domain.

Conclusions
We developed a new multiscale chemotactic motility model to investigate the behavior of ECPs in dynamic fluidic environments with spatially and temporally-varying gradients of two chemoattractants in response to ECP motility. We quantified the behavior of ECPs mimicking E. coli chemotaxis in fluidic environments containing the unmixed amino acid attractants N-methyl-L-aspartate and L-serine which function as strong chemical signals in the chemosensory system of E. coli. This was accomplished by formulating a novel dynamically coupled numerical model, MRC-CLB-ADT model, to simulate the motility of ECPs in an initially quiescent fluid with spatially and temporally-evolving chemoattractant concentrations. The MRC-CLB-ADT is capable of simulating the motility of ECPs by accommodating (i) spatial and temporal distributions of two distinct, non-interacting chemoattractants, (ii) effects of ECPs motion on the spatial and temporal distributions of the chemoattractants, and (iii) interactions of ECPs and the surrounding fluid environment.
The results and analysis of a variety of simulation set-ups allowed quantitative assessment for how chemoattractant distribution, particle-fluid dynamics, and particle-chemoattractant concentration field interactions affect the chemosensing properties of the ECPs. Results from "Imposed" simulations supported previous findings indicating that the ECP behavior is governed by receptor sensitivity rather than absolute attractant concentration [28,54]. Similarly, more recent biochemical work [55] has also substantiated the long-held notion that the networked architecture of chemosensory receptor arrays in E. coli plays a vital role in the cell's robust response behavior. Results of simulations that incorporate CLB-ADT emphasize the importance of advective-diffusive transport phenomena in modeling ECP trajectories in fluidic environments. Fluid effects on ECP chemosensing are significant both in those environments with pre-established chemoattractant gradients and in environments with distal point source chemoattractants. An example is shown in Appendix A.3: Figures A5  and A6 represent a more complex environment with solid obstacles and one attractant. ECPs will behave differently depending on whether or not chemoattractant is continually introduced into the environment. Further development and refinement of the model will be valuable for the exploration of other aspects of fluid-based ECP motility, such as the effects of fluid viscosity, non-Newtonian fluids, MCP stoichiometry in chemosensory arrays, chemoeffectors with varying affinities for MCPs, including repellent molecules, and environments comprising varying ECP population densities.

Acknowledgments:
The authors thank to Trinity University for awarding D. K. the Mach fellowship and providing computational resources for this project. The authors are grateful to Timothy R. Ginn of Washington State University, Pullman, WA for pointing us to the benchmark problem by [56].

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A Appendix A.1. Numerical Validation of the ADT Model
Prior to coupling with the MRC and CLB models, the ADT model was validated using a two-dimensional analytic solution provided in [56]. The analytic solution calculates spatiotemporal evolution of an inert substrate after being instantaneously injected into a Coutte flow at the center, x c = (x c , y c ), of the flow domain. The concentration distributions across the flow domain are computed by where x is the Cartesian coordinates of the lattice node, t is the time, Υ (t) is the variance matrix, det is the determinant operator. The variance matrix, Υ (t), is defined as where D 1 and D 2 are the diffusion coefficients and α = u x (x 2 )/x 2 , in which u x is the x-component of the specified fluid velocity at the top boundary and x 2 is the vertical distance between the substrate release location and the top boundary. The benchmark problem is illustrated in Figure A1. The model parameters in the LBM simulation were left in non-dimensional lattice units to be consistent with the benchmark problem. The size of the flow domain was set to 100 × 100 (lattice node numbers run from 1 to 101 in both directions). The fluid across the flow domain was initially stationary (i.e., u = 0). Specified fluid velocity was imposed at the top and bottom boundaries, but in the opposite directions. The lateral boundaries were assumed to be periodic. This resulted in zero flow velocity at the center of the flow domain (at x c = (x c , y c )). In the analytic solution, a point-like unit substrate concentration was instantaneously injected into the flow domain at the center, x c . The diffusion coefficient, D, of the substrate was set to unity (D = 1). The D 1 and D 2 in the variance matrix in Equation (A2) were also set to unity for an isotropic and homogeneous flow domain. From the half length of the domain, x 2 = 50, and the specified velocity at the top boundary, u(x 2 ) =1, we compute α = 2 × 10 −2 . Hence, the typical timescale associated with the shear rate is τ s = α −1 = 50.
Instantaneous point injection of the substrate in the analytic model [56] was represented by a point-like initial distribution of the substrate described by C (x c , t = 0) = δ (x c ), for which C → ∞ for t = 0 in Equation (A1). We simulated instantaneous substrate injection as a point source at x c in the LBM by implementing C (x c , t = 0) = 1, which resulted in ∑ x C (x, t) = C (x c , t = 0) = 1 for t > 0.
A comparison of the spatial distributions of the substrate computed by the analytic solution and numerical simulations via the ADT model for t = τ s and t = 2τ s are shown in Figures A2 and A3. Figure A2 shows that analytically and numerically computed spatial distributions and temporal evolution of the substrate are the same for t = τ s and t = 2τ s , except that the peak concentration at x c appears to be higher from the ADT model solution than the analytic solution at t = 0.2τ s . This was further confirmed from the horizontal and vertical cross-sections of the concentration field with respect to the center of the flow domain in Figure A3. Use of C (x c , t = 0) = δ (x c ) in the analytic model and C (x c , t = 0) = 1 in the ADT model led to discrepancies in the peak concentration, C (x c , t > 0) at early times while matching the concentration profiles almost perfectly away from the point source. At later times, concentration profiles computed by the ADT model matched the analytically computed concentrations perfectly. The total mass was strictly conserved in every time-step in the ADT model simulation.   Figure A4 demonstrates ADT model simulation of advective diffusive transport of a substrate (chemoattractant) in a horizontal flow channel with and without an internal obstacle. A no-slip boundary condition was implemented on the channel walls and a periodic boundary condition was imposed at the inlet and outlet. A steady Poiseuille flow was established in the horizontal channel prior to start of ADT model simulations. The Péclet number, Pe = u ss W/D was set to 500 in these simulations, in which u s is the average steady fluid velocity at the inlet prior to release of the chemoattractant into a flow domain and W is the channel width. The dimensionless time, Ω, was rendered as Ω = tu s /L, where L is the channel length and Ω = 0 corresponds to the time at which chemoattractant was released into the steady flow field. τ = 0.501 in these simulations. The results show that the obstacle enhances chemoattractant diffusion, which in turn alters its temporal and spatial gradients, which would have significant effects on motility behavior of a chemosensory particle. In Figure A5, the chemoattractant concentration was initialized at the center of the domain at t = 0 with no additional chemoattractant releases for t > 0. Hence, the total chemoattractant mass remained constant throughout the simulation. On the other hand, for the simulation results shown in Figure A6, the chemoattractant concentration was initialized at the center of the domain at t = 0 and in each subsequent time-step, t > t, C was released from the center of the flow domain.
In both simulations, the chemosensory particle, shown as a black circle, was initially placed outside the concentric ring of obstacles while the chemoattractant was initialized inside the ring. In MRC-CLB-ADT model simulations, the chemosensory particle sensed and avoided inline obstacles while navigating towards the spatially-and temporarily-varying maximum concentration gradients inside the ring. Consideration of the effects of inline obstacles on the spatial and temporal distributions of the chemoattractant concentration field, also simultaneously altered by particle-fluid hydrodynamic forces, is the unique modeling capability of the MRC-CLB-ADT model. (e) (f) Figure A5. Snapshots of the chemoattractant concentration field (whose contour lines are shown by solid lines) and trajectory of a chemosensory particle (shown by a light green dashed line) at time-steps of 1000; 28, 000; 68, 000; 80, 000; 174, 000; and 204, 000 in (a) through (f), respectively. The total chemoattractant mass remained constant throughout this simulation.
(c) (d) (e) (f) Figure A6. Snapshots of the chemoattractant concentration field (whose contour lines are shown by solid lines) and trajectory of a chemosensory particle (shown by a light blue dashed line) at time-steps of 1000; 14, 000; 32, 000; 42, 000; 52, 000; and 62, 000 in (a) through (f), respectively. Because the chemoattractant was continuously added to the center of the domain, the total chemoattractant mass in the fluidic domain increased in time. As a result, the domain was saturated with the chemoattractant at later times, which led to the random walk in particle trajectory towards the end of the simulation.  Scaling factor for methylation [17] 0.0625 b Scaling factor for demethylation [17] 0.0714 k Z Rate constant [17] 30/[CheZ] tot µM −1 s −1 k Y Rate constant [17] 100 µM −1 s −1 k s Scaling coefficient [17] 0.45 µM