Model

Superparamagnetic nanoparticles are magnetized in the presence of magnetic fields, and thus a slowly alternating magnetic field results in the oscillation of the particles. In our wsc1 activation, we attached 250 nm Ni-coated magnetic nanoparticles (MNPs) to his-tagged wsc1, and applied an external magnetic field to make MNPs oscillate. The oscillation will generate mechanical force on wsc1, which could be detected and induce the activation of the cell wall integrity (CWI) pathway. To calculate the actual force generated on wsc1 by magnetic activation, we constructed a model to analyze the net force experienced by MNPs under magnetic fields with different field strengths.

Magnetic force

In a simplified model, the magnetic force on a single magnetic nanoparticle (MNP) depends on three factors:

Volume of the MNP

The difference of magnetic susceptibility between MNPs and the medium

The magnetic field gradient

The magnetic force exerted on a single MNP could be calculated as (Subramanian et al,2019):

Based on Maxwell's equation, the equation could be further simplified as:

The calculation of the magentic force depends on the magnetic susceptibility difference between MNPs and the medium. Here is the model; we take water as the medium as our MNPs are suspended in water.The magnetic susceptibility of water is -9^-6, and the magnetic susceptibility of the single MNP is calculated as 1.34375. The difference of magnetic susceptibility is calculated as:

The 3-dimensional magnetic field gradient, Bx / By / Bz, will be obtained from the magnetic field simulation.

Drag force

The 250 nm MNP we will use is suspended in water. When it moves across the magnetic field, it will experience a drag force caused by the medium (in this case, water).

We use a simple equation (Stoke's law) for the calculation of drag force experienced by MNPs in water:

Where:

𝜂 is the viscosity of water, 1,0006 mPa·s

𝑅ℎ is the radius of the MNP. In this case, we assume the hydrodynamic radius of the particle equals to the MNP radius.

Δ𝑣𝑚 is the velocity of an MNP when it is moving in water

In our initial testing on MNPs, we took a video of MNPs moving in water in a specific direction. A small N52 neodymium magnet was placed next to the slide. This way, we could calculate the travelling speed of MNPs in water when a magnetic field is present. Four particles were selected and their velocity were calculated: the maximum is 0.7 um/s and the minimum is 0.4 um/s.

Video shows MNPs moving leftwards in water under the external magnetic field.

Net force calculation

The net force experienced by a single MNP is the difference between the magnetic force and the drag force. It is calculated with a simple method:

where Fn is the net force, Fm is the magnetic force, and Fd is the drag force caused by the medium.

A graph is plotted to see the change of Fn, Fm, and Fd as the field gradient changes. The numerical data we used here are hypothesized.

The model informs us that both the magnetic force and drag force will increase when the field strength increases, and the net force is the actual force that will be applied to wsc1 by MNPs. In our fluid-test video, we noticed that apart from the general movement towards left, MNPs also displayed brownian motion-like micromovement. This makes the calculation of the speed difficult as it is hard to tell if MNPs are moving in a constant speed or accelerating. However, we deduce that at a specific field strength, magnetic force will have the same magnitude as the drag force, and MNPs will move at a constant velocity.

In future, we are planning to combine our Magnetic Force model with our Magnetic field model for the coil to calculate the actual force experienced by a single MNPs under different magnetic fields. The model will also be corelated with wsc1 model to demonstrate the activation level of CWI pathway under different magnetic fields.

In order to find the exact value of how much force the magnet can apply on the magnetic nanoparticles, we ran a FEMM simulation, and then extracted the magnetic field gradient from the simulated result.

The model was obtained by running fine element simulation on FEMM software. The model was a replica of the Cavendish coil we used to generate MF, with peak magnetic field strength at ~130mT.

The simulation shows that we can generate sufficient MF gradient to activate our cells with MNPs. When we are designing the coil or activating the Wsc1 pathway, we should put the cell colonies at the edge of the iron core, which can maximize the gradient and the force applied to the MNPs.

When we apply a high frequency oscillating magnetic field to our media containing superparamagnetic nanoparticles, they begin to heat. Over time this heat generated diffuses to the surrounding media. For our goal of spatially targeted heating, it was essential to understand how the rate of heat generation compared to the rate of heat diffusion. To maximise our spatial resolution, we needed to be able to heat our small region before the heat equilibrated with the surrounding media.

To begin, we modelled just the heat generation of the MNPs. We know the spatial distribution of field strength can be approximated as a gaussian with a given bandwidth. If we know the heating efficacy at this peak the heating rate at any other point is given by: \[ SAR(x) = \frac{B(x)}{B_{max}}^2 SAR_{MAX} \]

We planned to disperse these MNPs in agar to determine how the bulk media would heat. For a given SAR, the temperature increase is given by: \[ \Delta T = SAR * /c_p \]

For a given heat distribution, the diffusion equation is: \[ \frac{\partial u}{\partial t} - \alpha \nabla u = \beta \partial t \]

Where temperature, \(u ~ f(x, y, y)\) subject to some heat generation \(\beta ~ f(x, y)\) calculated above and \(\alpha\), the thermal diffusivity.

Analysis of this can only be achieved using numerical methods. The updated rule is via forwards difference for x position, i, y position j and timestep k is: \[ u_{i, j}^{k+1} = \gamma (u_{i+1, j}^{k} + u_{i-1, j}^{k} + u_{i, j+1}^{k} + u_{i, j-1}^{k} - 4u_{i, j}^{k}) + \beta \Delta t + u_{i, j}^{k}\quad \text{with} \quad \gamma = \frac{\alpha \Delta t}{\Delta x^2} \]
This equation is stable for \( \Delta t \le \frac{\Delta x^2}{4 \alpha} \)

Analysing this for a plate initially at room temperature with a central peak SAR of 0.63W/g yields the following result:

Assumptions

Field Strength distribution is a gaussian in the x direction, with no variation in the y direction.

The surroundings remain at room temperature

The concentration of MNPs is homogenous

No heat is lost through convection or radiation

Conclusion

From this model, we were able to estimate the heat distribution we expected across our heated plates. This informed us that with our coil design, a high resolution of spatial targeting was feasible. A Hill function (in the pathway models) applied onto the approximately Gaussian distribution of the heat model output would reduce the number of cells that would be activated, thus resulting in higher precision.

To determine how this may impact our downstream pathways, we used the values generated from this model and implemented them as the gradients of a temperature function in the pathway models. In the TcI (and TlpA) and HSF model (iteration 4), we modelled fluorescent protein concentration against time for different temperatures at different positions of the magnetic field.

A notebook containing a python implementation of this analysis can be found in our Software Report.

Thermal: E. Coli (TlpA and TcI)

For TlpA and TcI it was crucial for us to understand how our constructs expressed TlpA and TcI and how this expression varied under different temperatures. Additionally, it was important to understand how TlpA and TcI bound to their respective promoters and repressed transcription, and how different temperature affected binding affinity.
Both TlpA and TcI operate in very similar ways, as they are both dimerising transcriptional repressors. Therefore both pathways were modelled with the same equations, and only their parameters varied.

Model 1.1 - Temperature Dependent ODE for TcI and TlpA

In this model we introduced temperature dependence into our ODE model specifically to influence how the DNA was transcribed in an attempt to emulate the function of TcI and its impact on transcription. We incorporate this temperature dependence in the form of an adapted hill function that utilises temperature as a variable, such that at higher temperatures transcription was uninhibited and at lower temperatures it was inhibited. We used a simple logarithmic function to simulate how temperature might vary in our hardware system such that it gradually asymptotes towards 42℃.

Assumptions

Assumed that temperature varies logarithmically

Our model assumes that TcI is already in excess

Our parameters were based on standard E. coli transcription and translation numbers; even though we expect ours to be largely the same there may be some variations

We had to make estimations on what the maximum rate of transcription would be and at what temperature half the maximum binding of TcI/TlpA would be achieved

Conclusion

From this model we can see that it is possible for our constructs to start expressing the fluorescent protein in a reasonable amount of time (approximately 1 hour), and gave us a good indication of the minimum time needed to heat our samples in wet lab experiments to see a visible change. Taking into account the long maturation time of mRFP1 (1h), we then checked the samples both after an hour as well as several hours later.

Model 1.2 - Different temperatures in the ODE

This iteration of the model incorporated temperature gradient values from the heat modelling in linear temperature equations, which assumed no heat loss. This showed the difference in expression depending on the temperature, and thus the position within the magnetic field.

Conclusion

We achieved simple ODEs to represent the molecule pathways in the TlpA and TcI constructs. We also began to combine them with the heat model, in order to predict the level of fluorescence depending on the position in the magnetic field through a temperature function.

Our future iterations would have included more refined temperature models, which were non-linear and included diffusion of heat. This would help inform us where the area of highest activation on the petri dish is. This is important in pathways such as TcI, where the highest temperature may not result in the most fluorescence (see AIS-China 2023’s characterisation of part BBa_C0051). Therefore areas closest to the coil may not result in the most fluorescent output. The combination of these models would also inform us of the levels of pathway activation beyond our initial targeted area, thus feeding back to the DBTL cycle of hardware to increase precision of the spatial targeting.

Combination of the pathway and toxicity models would help inform the optimum concentration of MNPs that balance the amount of heat produced as well as population growth. Therefore, this would guide the concentration of MNP in agar plates made for wet lab experiments.

Model 2 - Temperature Dependency of TcI’s Binding Affinity

The aim of this model was to find the ideal temperature where we could achieve sufficient transcription of our fluorescent protein without subjecting our cells to high temperatures that might lead to death at prolonged periods of time. In this model we estimated the binding affinity TcI by normalising our data to fluorescence at room temperature (maximum binding/repression) and fluorescence at 45℃ (minimum binding/repression). We defined the critical temperature (or optimal temperature in our case) as the point where the curvature of our sigmoid changes. It is important to highlight that the data used in this model is purely simulated, but would have been real if time allowed.

Conclusion

With more time and data, we would have improved the parameters in this model, and fed it into the expression model 1.2 and heat model to guide hardware designs. Additionally, we would have compared the two different model approaches for TcI (and TlpA), and determined which model was more suitable for the pathway depending on the wet lab data obtained.

To find more details on the models and parameters used:

Thermal: S. cerevisiae (HSF)

Heat shock factor (HSF) is a transcription factor native to yeast that upregulates transcription of various heat shock proteins (HSPs) when temperature increases (Mühlhofer et al., 2019). HSF trimerises, and during heat shock, transcription is increased through two different models: chaperone titration model and phosphorylation model (Zheng et al., 2016). In our design of the HSF regulated construct, the HSP26 promoter which is regulated by HSF is upstream of a fluorescent protein (FP). This model aimed to show the change in FP over time depending on the temperature.

Transcription: \[\frac{d[mRNA]}{dt}=k_{1}[DNA]\frac{T^n}{K+T^n}-λ_m[mRNA] \]

Translation of dark FP: \[\frac{d[P_d]}{dt}=k_{2}[mRNA]-λ_{pd}[P_d] \]

Maturation of dark FP: \[ \frac{d[P_m]}{dt}=k_3[P_d]-λ_{pm}[P_m] \]

In iteration 4 of the model, the concentration of mRNA, dark protein and mature protein is modelled as an ODE against time at varying temperatures. These temperatures are gradient values informed by the heat modelling under the assumption that there is no heat loss (for other iterations and parameter assumptions, see pdf below).

Conclusion

We wanted to use this model to inform the team of the difference in rate of fluorescence as well as minimum incubation time depending on the position of the cells within the magnetic field through the temperature function. Unfortunately, we did not have a working construct for this pathway, therefore it was not used to inform our wet lab experiments.

However, future iterations would have involved:

More accurate parameters informed by wet lab data

A non-linear temperature function that considered the dissipation of heat

An improved ODE such that a constant temperature of 25℃ would not result in an increase in translation and transcription

Combined heat, pathway and toxicity modelling to guide the wet lab team on the optimum temperature and MNP concentration for maximum growth and fluorescence

We would then use this model to feedback into the hardware design and try to maximise the difference in temperature gradient across a plate. This would guide us on how to maximise the difference in fluorescence between the targeted and non-targeted regions.

See below for more details on iterations of the model, parameters and assumptions:

Mechanical (Wsc1)

Wsc1 is a cell membrane protein in S. cerevisiae which acts as a sensor for the cell wall integrity pathway, activating the Rom2 GEF and thus the downstream MAPK signalling cascade to regulate the expression of at least 25 genes implicated in cell wall biogenesis (Philip and Levin, 2000). This is done partially through the Rlm1 transcription factor, which we have used as our reporter in the wet lab designs.

The cell wall integrity pathway is complex to model due to the numerous steps (figure 6), most of which have not had their constants characterised. This model therefore treats the activation of Rlm1 by Wsc1 as a single step, when in reality it involves a complex pathway including a MAP Kinase cascade. This abstraction relies on the assumption that only one of the steps in the cascade is rate-limiting. If multiple steps operate on a similar time scale, there will be a lag in the response not captured in this model. This could be modelled in the future using delay differential equations.

In our model, Rlm1 transcription is dependent on the proportion of activated Wsc1. Wsc1 activation itself is modelled as a function of “magnetic flux”. This measure of magnetic flux is entirely arbitrary, and an interesting future step would be to characterise the relationship between field strength/gradient.

Conclusion

Turning the activation and protein expression ODEs into a population level model would allow for the integration of the effects of toxicity. This would enable us to identify what concentration of magnetic nanoparticles is high enough to achieve strong induction, while low enough to maintain population levels.

Additionally, this Wsc1 pathway model could be combined with the MNP drag force model to predict which areas of the petri dish would be activated the most, thus feeding back into the hardware design to increase the precision of targeting.

To find more details on the models and parameters used:

Magnetic (EPG)

EPG is a novel protein that exhibits sensitivity to magnetic fields, although the mechanism behind this property is not yet well-understood. The kinetic model of EPG remains largely unexplored, and the optimal method for activating EPG in our setup is currently unclear.

As part of our research, we are designing new fusion proteins, and we aim to use in-silico computations to identify the optimal effector protein and appropriate linker sequences. To this end, we used AlphaFold3 to simulate the structures of various constructs, including EPG with and without its signal peptide, EPG-LacI fusion proteins with different linkers, and fusion proteins involving EPG-TEV and EPG-Nanoluciferase.

AlphaFold's capacity to predict different conformations of proteins based on their energy levels was utilized to simulate the potential states of EPG and its fusion proteins. Since AlphaFold predicts structures that correspond to favorable energy states, this enabled us to explore how the protein might fold under various conditions. However, we understand that AlphaFold does not provide molecular dynamics simulation of the exact protein structure. It can only provide a statistically sampling of the protein structure. However, for an initial demonstrative model, it is accurate enough for us to know how the EPG can behave.

Our initial construct was an EPG-LacI fusion protein, where EPG was inserted into the second loop of LacI to potentially regulate LacI activity via magnetic field sensing. However, the long signal peptide of EPG proved to be problematic. The presence of this peptide and the linker sequence caused spatial interference with the DNA-binding domain, leading to a non-functional protein. Additionally, since the site of insertion was meant to act as a hinge, the inherent hinge-like nature of EPG limited its ability to sense magnetic fields effectively. Consequently, the EPG-LacI system was deemed unsuitable and was discontinued.

To further understand the folding and hinge-like activity of EPG, we aligned the structures of EPG with and without its signal peptide. The alignment revealed that a specific region within EPG facilitates the movement of the signal peptide, possibly enabling signal relay. Furthermore, we found that the core Ly6 domain of EPG is highly conserved and can exist in two distinct states, potentially corresponding to its active and inactive conformations. This suggests that magnetic fields might induce specific conformational changes in the Ly6 domain, toggling between folded and unfolded states.

Next, we evaluated EPG-TEV fusion proteins with different linkers. When a rigid 5 amino acid (aa) linker (PAPAP) was used, the structure displayed distorted signal peptide folding in its closed conformation. The rigidity of the linker restricted conformational space and forced the protein into a specific orientation, hindering activation. Additionally, the short signal peptide prevented the two split TEV fragments from coming into contact. A simulation using a longer, 10aa rigid linker revealed a similar issue, though the increased linker length reduced the strain within the protein and potentially lowered the activation barrier. However, there was concern that the longer linker might result in constant activation of EPG, leading us to favor shorter linkers, specifically the flexible 5aa (F5) and rigid 5aa (R5) linkers.

For the EPG-TEV construct with an F5 (SGGGS) flexible linker, we observed that the lack of spatial strain allowed for seamless binding of the TEV fragments. However, we suspect this construct might produce false negatives, as the absence of a spatial barrier could lead to TEV activation even in the absence of a magnetic field. Therefore, while the flexible linker enables interaction between the TEV fragments, it may not sufficiently regulate activation based on magnetic field presence.

Since both constructs show promising potential, we decided to prioritize implementing these two in our engineering cycle. Additional constructs will be produced and screened as time permits. The final construct will aim to balance the tension between structural strain and the ease of fragment joining. Achieving this balance will require further modeling or large-scale screening of linker sequences.

---

Model on Expression Levels

The expression level model serves as a preliminary guide to estimate the amount of protein produced within the cells, assisting in the selection of an appropriate promoter for expressing EPG-Nanoluc and EPG-TEV. Since these recombinant proteins are expressed in a heterogeneous host, it is crucial to prevent the formation of inclusion bodies. Additionally, achieving an optimal sensitivity and signal-to-noise ratio requires balancing the expression levels of the reporter and EPG-TEV constructs.

The model indicates a linear relationship between protein concentration and promoter activity. Therefore, instead of using the strongest constitutive promoter, we opted for the second strongest, J23100, to minimize cellular burden and reduce inclusion body formation.

This model also informs the timing of protein expression in our protocol. We found that DH5α cells have a doubling time of approximately 25–30 minutes and reach stationary phase in ~18 hours, while BL21(DE3) cells grow faster with a doubling time of ~20 minutes, reaching stationary phase in ~12 hours. Since our system is expressed in BL21(DE3) cells, and protein expression begins around 20 minutes after reaching steady state, we incubate the culture for 24 hours at 30°C to achieve the highest protein concentration.

---

Model on EPG-Nanoluc Activity Under Magnetic Activation

The second aspect of the model predicts the activity of EPG-Nanoluc when exposed to magnetic fields. We observed a decrease in Nanoluc activity, contrary to what has been previously reported in the literature. To better understand this discrepancy, we constructed a Delay Differential Equation (DDE) model, incorporating an intermediate state that accounts for the possibility of EPG forming dimers or multimers. This cooperative behavior was modeled using a Hill function, while the DDE framework accounted for the lag phase between magnetic activation and the observed change in activity.

When overlaying the activation curve onto a background activation pattern, the model suggests that during periods of negative EPG production—when cellular activity is reduced and EPG is being degraded—there is a decrease in EPG levels upon removal from the incubator. This degradation explains the observed dip in activity, as EPG breakdown outpaces its synthesis.

The model provides a feasible explanation for the reduction in activity seen after magnetic activation, as it indicates continuous degradation of EPG over time could lead to the dip on the background activation. However, due to current limitations in characterization techniques, we can only measure the decline in activity post-activation.

In future work, improved techniques capable of accurately measuring luminescence under magnetic field application would enable detailed characterization of EPG’s response across various magnetic field intensities and activation durations. This, in turn, would allow for further optimization of the model’s parameters and schematic design.

We aimed to construct a model of the relationship between magnetic nanoparticle toxicity and the values of growth parameters. While we have created a pipeline that can estimate growth parameters from growth curves under a range of magnetic nanoparticle concentrations, we were unable to obtain sufficient wet lab data to populate a model. Ideally, this would have been integrated into pathway ODEs at a population level to see how to balance induction with health of cells.