- Modelling -

Equations

In our first model, we utilise a system of ordinary differential equations representing the bacterium's gene regulatory network, in order to inform parameters for optimal uptake and cell stress.

We first defined a function of time representing the production of LutH, composing two Hill equations—the first representing the natural LutH production as exists in the wild-type, which we assume to be inhibited by lanthanide concentration, and the second representing our engineered LutH production, induced by IPTG:

We then defined a function representing the production of tLanM. We define P_H(t) and P_M(t) to use different IPTG concentrations, as ideally the promoters would use different inducers. The function composes a single Hill equation, representing the IPTG-induced tLanM production:

Using these two production functions, we form a function representing the stress in the system. Stress is increased by lanthanide ion concentration and greater production of LutH and tLanM. As the lanthanide concentration approaches L_d, stress approaches infinity:

We then use this function in χ(t), representing the proportion to which he rate of LutH and LanM production are scaled:

We then represent the formation of tLanM-Ln complexes by a proportion of the total tLanM and Ln present in the system:

Rate of change of LutH is then represented by a differential equation, whereby production of LutH is scaled down by stress, while degraded LutH is removed from the system at a constant rate:

Rate of change of tLanM is implemented similarly, with tLanM production being scaled down by stress, then decayed tLanM being removed. We also subtract the tLanM in tLanM-Ln complexes.

Rate of change of lanthanides in the cell is represented in 3 parts, composing the subtraction of decayed lanthanides and those sequestered by the tLanM, with a term representing uptake from the environment. We use the Michaelis-Menten enzymatic kinetics function as we assume that the mechanics will be largely similar, scaled by the amount of LutH in the system:

Population Model

We additionally wanted a model to simulate the entire population across a larger space, in order to be able to visualise our system's behaviour on a much grander scale, potentially even to the scale of an industrial implementation, and also to inform us as to the parameters of our experimental design. To do this we used a system of partial differential equations.

Parameters

With consistent units (ommitted parameters have been determined indirectly by fitting the model with the wet lab data):

Hyperparameters

• Grid size. dx and dy are expressed in µm.
• Timestep size. Expressed in seconds.
• Lanthanides lag.

Assumptions

Grid and space assumptions

• The space is represented on a 2D square grid with uniform grid spacing, which assumes that bacterial motion, nutrient diffusion, and chemotaxis are isotropic and homogeneous across space.

Bacterial behaviour

• The bacterial population is modelled using a concentration field u, and it is assumed that this concentration follows diffusion and chemotactic behaviour.
• Bacteria are attracted to a chemotactic agent and nutrients via gradients in chemoattractant (v) and nutrient (f) concentrations.
• Bacteria undergo logistic growth based on the nutrient availablility and their maximum growth rate (r) up to a saturation point.
• The bacteria die at a constant rate (kd) in the absence of food, and they can take up lanthanides (l), which activates chemotaxis.

Chemotaxis sensitivity

• The chemotactic sensitivity to nutrients and the chemoattractant is a decreasing function of the cell concentration to account for overcrowding.

Lanthanides behaviour

• Lanthanides diffuse through the system with a constant diffusion coefficient DI and are taken up by bacteria with a rate proportional to the bacterial concentration (omega * u * l)
• The uptake of lanthanides has a delayed effect on chemotaxis, modelled by a lag period (lag) in time before the chemotaxis behaviour is expressed.
• The average uptake per cell (l_avg= 50) is considered constant.

Nutrient behaviour

• Nutrient diffusion occurs with a constant diffusion coefficient Df, and it is consumed by bacteria at a rate proportional to the bacterial concentration.
• The nutrient consumption follows Michaelis-Menten kinetics, where the consumption rate saturates at higher nutrient concentrations (f/(h + f)) with half-saturation constant h).

Chemoattractant dynamics

• The chemoattractant concentration v diffuses according to a constant diffusion coefficient Dv without additional production or decay terms.
• The chemoattractant is initially set up with a localised peak, simulating a source of attraction in the bacterial environment.

Cell activation

• A subpopulation of active cells (a) is considered, which responds both to chemoattractant and nutrient gradients, and undergoes diffusion and death similar to the bacterial population.
• Active cells diffuse like the bacterial population but have a distinct behaviour influenced by chemoattractant uptake.

Time-stepping

• The simulation uses a fixed time-step size for all processes, assuming that the timescales of diffusion, chemotaxis, bacterial growth, and lanthanide uptake are comparable.

Boundary conditions

• Cell concentration (u and a) are constrained to have zero gradients on the boundaries (no-flux boundary conditions), meaning that no material flows in or out of the boundaries.

Initial conditions

• Randomised initial conditions for bacterial (u), nutrient (f), and lanthanide (l) concentrations, with a Gaussian distribution for chemoattractant (v), assume some initial variability or noise in the system.

Limitations

Simplified spatial representation

• 2D grid assumption: the model assumes a 2D space, while bacterial motion and chemotaxis typically occur in three dimensions. This could lead to inaccuracies when extrapolating results to real-world 3D environments.
• Uniform grid spacing: the grid spacing is uniform, which may not capture finer spatial variations in concentration gradients, especially in regions with steep changes.

Constant diffusion coefficients

• Fixed diffusion coefficients: the diffusion coefficients for bacterial concentration (Du), chemoattractant (Dv), nutrients (Df), and lanthanides (Dl) are assumed to be constant. In reality, diffusion rates may vary with environmental factors such as temperature, viscosity, and the presence of obstacles or other microorganisms.
• The model assumes isotropic diffusion (same in all directions), but in real systems, bacterial motion could be influenced by external forces or environmental heterogeneity, leading to anisotropic behaviour.

Over-simplified chemotactic response

• Constant sensitivity: the chemotactic sensitivity parameters for nutrients account only for overcrowding effect. In real bacterial systems, the sensitivity to chemoattractant and nutrients may change over time or depend on the bacterial state, adaptation, or saturation effects.
• No receptor saturation: the model assumes that the chemotactic response scales linearly with the chemoattractant gradient, without considering potential receptor saturation or desensitization at high chemoattractant concentrations.

Nutrient consumption and growth

• Simplified nutrient kinetics: the nutrient consumption follows a Michaelis-Menten-like relationship, which may not fully capture the complexity of bacterial metabolism, including the possibility of cooperative uptake or the presence of multiple nutrient sources.
• Constant growth and death rates: the maximum growth rate (r) and death rate (kd) are constant, while in real bacterial systems, growth rates can vary depending on nutrient availability, population density, environmental stress, and quorum sensing.

Lanthanide uptake

• Lag period: the model assumes a fixed lag period (lag) between lanthanide uptake and chemotaxis adjustment, but in reality, this delay could vary dynamically depending on the concentration of lanthanides or the physiological state of the bacteria.
• No lanthanide-dependent feedback: the model does not account, apart from chemotaxis, for any feedback mechanisms where lanthanide uptake could influence other aspects of bacterial behaviour, such as metabolism or motility, beyond the fixed uptake rate.
• Constant lanthanides concentration: every cell is modelled uptaking a constant number of lanthanides that is released after the cell death.

Boundary conditions

No-flux boundaries: the zero-gradient (no-flux) boundary conditions prevent any flow of material in or out of the system, which may not accurately represent real-world systems where bacteria or chemicals could diffuse beyond the boundaries of the modelled region.

Absence of environmental factors

• No environmental variations: the model does not account for external environmental factors like fluid flow, temperature gradients, or mechanical forces, which can influence bacterial motility and chemotaxis in natural or industrial settings.
• No interaction with other species: the model focuses only on a single bacterial population, but real environments often involve interactions with other microorganisms, bacteriophages, or physical structures, which could alter chemotactic behaviour.

Simplified active cell dynamics

• No cell-state transitions: the model assumes two distinct populations—active cells (a) and non-active cells (u), but it doesn't account for transitions between these states based on environmental conditions or intracellular signalling.
• Fixed yield coefficient: the yield coefficient (Yc) for the conversion of nutrients into biomass is constant, which ignores variations in metabolic efficiency based on nutrient quality or bacterial energy state.

Absence of detailed biophysical processes

• No intracellular dynamics: the model focuses on population-level behaviour (concentration fields) but ignores detailed intracellular processes like gene expression, receptor dynamics, and energy metabolism, which play key roles in bacterial chemotaxis and growth.
• No cell-cell interactions: the model does not account for direct interactions between individual bacteria, such as collision, quorum sensing, or biofilm formation, which can significantly alter chemotactic and growth dynamics in dense populations.

Time-stepping and numerical stability

Fixed time step: a fixed time step is used for solving the partial differential equations (PDEs). While convenient for simulations, it may lead to numerical instability if the time step is too large relative to the rate of diffusion or chemotactic response, especially in regions with steep gradients.

Mathematics

Variables

Equations

---

1. Bacterial Population:

   $\frac{\partial u}{\partial t}=g(u)+\nabla ·(D_u\nabla u-\lambda (u)\nabla f)-\tau (u,l(x,y,t-\mathrm{lag}))$

   Where:

   • g(u) is the growth function.
   • D_u is the diffusion coefficient for bacterial cells.
   • λ(u) is the chemotactic sensitivity to nutrients.
   • τ(u,l(x,y,t-lag)) represents the transition to the active state when chemotaxis is enabled.

---

2. Chemoattractant:

   $\frac{\partial v}{\partial t}=k(u,v)+\nabla ·(D_v\nabla v)$

   Where:

   k(u,v) models the kinetics of the chemoattractant
   D_v is the diffusion coefficient of the chemoattractant

---

3. Nutrients:

   $\frac{\partial f}{\partial t}=\nabla ·(D_f\nabla f)-\frac{1}{Yc}r\left(\frac{f}{h+f}\right)(u+a)$

   Where:

   • D_f is the diffusion coefficient for nutrients.
   • Yc is the yield coefficient for cell growth.
   • r is the maximum growth rate

---

4. Lanthanides concentration:

   $\frac{\partial l}{\partial t}=\nabla ·(D_l\nabla l)-\omega ul+k_d{l}_{\mathrm{avg}}$

   Where:

   • D_l is the diffusion coefficient for lanthanides.
   • ω is the rate at which cells uptake lanthanides
   • k_d is the death rate of cells in the absence of nutrients
   • l_avg is the average amount of lanthanides per cell

---

5. Active cells

   $\frac{\partial a}{\partial t}=\nabla ·(D_u\nabla a-\lambda \nabla f-\chi (v)\nabla v)+\tau (u,l(x,y,t-\mathrm{lag}\left)\right)-k_{d}a$

   Where:

   • χ(v) is the chemotactic sensitivity dependent on the chemoattractant concentration.

Description

The model provided for engineered chemotaxis of Methylobacterium extorquens AM1 after lanthanide uptake is fundamentally an extension of the classic Keller-Segel model for chemotaxis. The Keller-Segel model describes how cells (such as bacteria) move in response to a chemical signal (the chemoattractant) in their environment.

The Keller-Segel model is composed of two primary partial differential equations:

• Cell density u(x,t): describes how the concentration of cells changes over time due to diffusion and chemotaxis.



$\frac{\partial u}{\partial t}=\nabla ·(D_u\nabla u-\chi (u,v)\nabla v)+\mathrm{growth\; term}$



• Chemoattractant concentration v(x,t): describes the concentration of the chemoattractant that cells are attracted to, which also diffuses through space and may be produced or degraded.



$\frac{\partial v}{\partial t}=\nabla ·(D_v\nabla v)+\mathrm{production}\bigvee \mathrm{decay\; term}$

Where:

• D_u is the diffusion coefficient of the cells.
• χ(u,v) is the chemotactic sensitivity, a function that models the strength of chemotaxis.
• Δu represents the diffusion of cells.
• ∇⋅(χ(u,v)∇v) is chemotaxis, or movement in response to the chemoattractant gradient.
• D_v is the diffusion coefficient of the chemoattractant.

Key differences

While the Keller-Segel model focusses on simple bacterial chemotaxis to a chemical signal, this model introduces additional complexity:

• Chemotaxis towards nutrients: we include chemotaxis not only towards the chemoattractant, but also to nutrients, coupling bacterial movement with nutrient availability.
• Lanthanide uptake: the absorption and release of lanthanides influences cell behaviour, introducing a delay (lag) before chemotaxis can occur after lanthanide uptake.
• Active cells: our model distinguishes between the general bacterial cells u, and actively chemotactic cells a, with the latter only exhibiting chemotactic behaviour after lanthanide uptake.

Simulation

In order to solve our PDE system, we used the Python library FiPy, a finite volume PDE solver developed by NIST.

We first initialised our parameters, initial conditions, and boundary conditions:

# Define the grid
nx = ny = 120
dx = dy = 100 
mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
np.random.seed(70)

# Define the parameters
Du = 1000      # Diffusion coefficient for u 200 to 1900 µm^2/s 
Dv = 600       # Diffusion coefficient for v 600 µm^2/s
chi = 1000     # Chemotactic sensitivity 10e3 to 10e5 µm^2/s
alpha = 0.1    # Production rate of v
beta = 0.05    # Decay rate of v
r = 0.00005    # Maximum specific grow rate (3.89 to 6.94)x10^-5 per second
c_max = 1      # Maximum concentration of cells
Df = 1200      # Diffision term for nutrients 600 µm^2/s
h = 1          # Half saturation constant 0.2 to 2.5 mM 
theta = 1000   # Chemotactic sensitivity for food
Yc = 0.44      # Yield coefficient gram of cells per gram of food 
kd = 0.0005    # Rate at which cells die in the absence of food (5.08±0.55)x10ˆ4 per second per 1000 cells
omega = 0.01   # Uptake rate of lenthanides scaling factor
Dl = 600       # Diffusion coefficient lenthanides 
lag = 5        # Lag between lenthanides uptake and chemotaxis towards chemoattractant
lenc = []      # Lanthanides concentration
max_v = 50     # Maximum chemiotractor concentration in mMol
max_f = 50     # Maximum initial concentration of nutrients in mMol
lavg = 50      # Average uptake of lanthanides per cell 50 mMol

chem_centre = 0 # Centre of the chemoattractant 

# Define the variables
u = CellVariable(name="u", mesh=mesh, value=1.0)       # Bacterial Concentration 
v = CellVariable(name="v", mesh=mesh, value=0.0)       # Chemoattractant Concentration 
f = CellVariable(name="f", mesh=mesh, value=0.0)       # Nutrients Concentration 
l = CellVariable(name="l", mesh=mesh, value=0.0)       # Lanthanides Concentration
a = CellVariable(name="a", mesh=mesh, value=0.0)       # Concentration of active cells 
diff = CellVariable(name="diff", mesh=mesh, value=0.0) # Difference between Lanthanides concentration

# Initial conditions (localized peak in the center)
x, y = mesh.cellCenters
x_centre, y_centre = nx/2, ny/2
r_squared = (x-x_centre)**2 + (y-y_centre)**2                  

# Set the initial conditions 
u.setValue(0.05 * np.random.rand(nx*ny)) # Initial bacteria concentration 0.01 - 0.1 bacteria/µm^2

v.setValue(max_v * np.exp(-((x - nx*dx/2 + chem_centre)**2 + (y - ny*dy/2 + chem_centre)**2) / (2*(2000**2))))

f.setValue(50)

# f.setValue( 0.5 * max_f * np.ones(nx*ny) + 0.5 * max_f * np.random.rand(nx*ny))
l.setValue(0.3 * np.random.rand(nx*ny))
lenc.append(l.value.copy())

a.setValue(0.0) # Initialize active cells to 0


# Set boundary conditions

# Impose low flux on the left and right boundaries (Neumann BC)
u.faceGrad.constrain(0, mesh.exteriorFaces)
a.faceGrad.constrain(0, mesh.exteriorFaces)
f.faceGrad.constrain(0, mesh.exteriorFaces)
l.faceGrad.constrain(0, mesh.exteriorFaces)

u.constrain(0.01, mesh.exteriorFaces)
a.constrain(0, mesh.exteriorFaces)
l.constrain(0, mesh.exteriorFaces)
	

Then defined our equations:

# Define the PDEs
eq_u = (TransientTerm(var=u) == 
        -((theta * (1/(1+(1+u+a)**2)) * u * f.grad[0]).grad[0]      # Attraction to nutrients
          + (theta * ((1/(1+(1+u+a)**2))) * u * f.grad[1]).grad[1]) #
        + ((Du * u.grad[0]).grad[0] + (Du * u.grad[1]).grad[1])     # Cell diffusion term
        + (r * (f/(h + f)) * (u+a) - kd/1000 * u)                   # Cell growth term 
        + diff * u/20)                                              # Cell Activation term
        

eq_v = (TransientTerm(var=v) == ((Dv * v.grad[0]).grad[0] + (Dv * v.grad[1]).grad[1])) # Chemoattractant Diffusion term 

eq_f = (TransientTerm(var=f) == 
        ((Df * f.grad[0]).grad[0] + (Df * f.grad[1]).grad[1]) # Nutrients Diffusion term 
        - 100/Yc * r * (f/(h + f)) * (u+a))                   # Nutrients consumption


eq_l = (TransientTerm(var=l) == 
        ((Dl * l.grad[0]).grad[0] + (Dl * l.grad[1]).grad[1]) # Lanthanides diffusion term 
        - omega * u * l                                       # Lanthanide absorption by cells 
        + kd/1000 * a * lavg)                                 # Lanthanides release by cells 

eq_a = (TransientTerm(var=a) == 
        -(((chi * (((max_v-v)/max_v)/(1+(1+u+a)**2)) ) * a * v.grad[0]).grad[0]    # Chemotaxis
          + ((chi * (((max_v-v)/max_v)/(1+(1+u+a)**2)) ) * a * v.grad[1]).grad[1]) # 
        -((theta * ((1/(1+(1+u+a)**2))) * a * f.grad[0]).grad[0]                   # Attraction to nutrients
          + (theta * (1/(1+(1+u+a)**2)) * a * f.grad[1]).grad[1])                  #
        + ((Du * a.grad[0]).grad[0] + (Du * a.grad[1]).grad[1])                    # Cell diffusion term
        - diff * u/20                                                              # Cell activation term
        - kd/1000 * a)                                                             # Cell death
        
	

And finally solved and visualised the system:

# Visualisation
for step in tqdm(range(steps), desc="Time-stepping"):
    
    vconvection_x = (chi * (((max_v-v)/max_v)/(1+(1+u+a)**2)) ) * a * v.grad[0]
    vconvection_y = (chi * (((max_v-v)/max_v)/(1+(1+u+a)**2))) * a * v.grad[1]


    f_grad = f.grad
    fconvection_x = theta * (1/(1+(u+a)**2)) * a * f_grad[0]
    fconvection_y = theta * (1/(1+(u+a)**2)) * a * f_grad[1]

    # Compute the divergence of the convection term
    div_convection =  - (fconvection_x.grad[0] + fconvection_y.grad[1])


    if len(lanc) >= lag:
        diff.setValue(lanc[-1] - lanc[-lag])
        diff.setValue(np.minimum(0,diff.value))


    # Updating Lanthanides Concentration 
    lanc.append(l.value.copy())
    if len(lanc) > lag:
        lanc.pop(0)
   
    if step % print_step == 0:
        with h5py.File(hdf5_file, "a") as r:
            r["cell"][int(step/print_step)] = u.value.reshape(nx, ny)  
            r["Ln"][int(step/print_step)] = l.value.reshape(nx, ny)  
            r["food"][int(step/print_step)] = f.value.reshape(nx, ny)  
            #r["Divergence"][step/print_step] = div_convection.value.reshape(nx, ny)  
            r["active cells"][int(step/print_step)] = a.value.reshape(nx, ny)
            #r["chem"][step/print_step] = v.value.reshape(nx, ny)  

    # Solve the PDEs
    eq_u.solve(dt=timeStepDuration)
    eq_v.solve(dt=timeStepDuration)
    eq_f.solve(dt=timeStepDuration)
    eq_l.solve(dt=timeStepDuration)
    eq_a.solve(dt=timeStepDuration)

    # Add Gaussian noise to variables
    u.setValue(u.value + noise_level_u * np.random.normal(0, 1, u.value.shape))
    v.setValue(v.value + noise_level_v * np.random.normal(0, 1, v.value.shape))
    f.setValue(f.value + noise_level_f * np.random.normal(0, 1, f.value.shape))
    l.setValue(l.value + noise_level_l * np.random.normal(0, 1, l.value.shape))
    
    u.setValue(np.maximum(u.value, 0))
    v.setValue(np.maximum(v.value, 0))
    f.setValue(np.maximum(f.value, 0))
    l.setValue(np.maximum(l.value, 0))
    a.setValue(np.maximum(a.value, 0))
    f.setValue(np.minimum(f.value, max_f))
        
print("Frames saved.")

	

Analysis

Running the simulation with the chemoattractant emanating from the centre of the plane, we produced the following results:

As time passes, the environmental lanthanide concentration (top right) decreases as the bacteria uptake the lanthanides, while the cells move towards the centre, bringing the lanthanides with them. This then also represents the net movement of lanthanides towards the centre.

Running the simulation with the chemoattractant emanating from the left side of the plane, we produced the following results:

This configuration was shown to be comparably effective, showing the cells moving towards the left side.

The population model significantly contributed to our experimental design and wet lab procedures by enabling predictions of assay outcomes under various conditions. It firstly allowed us to anticipate the behavior of the bacterial population in response to different gradients, helping us to forecast the results of future assays. Secondly, the model suggested that a central inoculation point for the chemoattractant would be more efficient for maximising bacterial chemotaxis. This insight directly informed the design of our experiments, ensuring more effective use of resources and time in the wet lab. Lastly, the model's ability to simulate different diffusion coefficients allowed us to experiment with varying media densities. This helped us select the optimal media concentration, improving bacterial motility and chemotactic response, leading to more reliable experimental results. Overall, the model was a powerful tool for both guiding experimental decisions and enhancing the accuracy of the wet lab's experimental outcomes.