Model

Nominated for best model, our hardware and biological systems are shown here is the documentation

The integrated modelling and simulation section of iGEM Sheffield 2023 describes the development of a computational representation of the function of the PARSE genetic system in Escherichia coli. Additionally, this system can be applied for the maintenance and control of bacterial cultures for experimental purposes.

The first section describes the general structure of the PARSE system, including the plasmids, genes and sequences that compose them, as well as the required information to ensure reliable monitoring of gene expression. The involved parameters and interactions that direct gene function are described in mathematical terms. Our experimental results are compared with the behaviour of the system under optimal conditions, and discrepancies are discussed.

Then, the detailed reaction mechanisms through which the PARSE system regulates growth rate by its selective expression are explored and the requirements and considerations for its potential application for experimental tasks are listed.

Finally, the design, features and operation of a turbidostat developed specifically to work with optimal performance in combination with the PARSE system to grow and monitor bacterial cultures in a controlled, precise and affordable manner, is introduced.

Biological Modelling

Hardware Modelling

Hardware Modelling

Introduction

The hardware modelling of the project was created to simulate and aid in specifying and validating the system requirements for a multiplexed turbidostat bioreactor. This system is intended for growth modulation experiments across various scales, offering a more streamlined and cost-effective alternative to intricate and expensive microplate spectrophotometers. This type of bioreactor makes it possible to automatically build an inducer vs OD standard curve, and to obtain several similar curves for growth characterisation purposes in parallel. With its several cells that may contain different microbial cultures, multiple cultivation experiments under varying conditions can be run simultaneously and analysed in real-time while maintaining turbidity at a constant level.

This section subdivides and offers representations of the fundamental controls for each subsystem, mirroring the individual control systems within the physical bioreactor. This aspect of the modelling acts as a link between the growth equations described in the biological modelling and the system hardware design. By creating the Simulink model, the hardware team is aided in understanding how, for example, the temperature sensor readings and the fluid system interact. Additionally, the gene circuits built in SimBiology dictate relevant parameters for the bioreactor to measure and detail equations and parameters. This informed construction requirements in the bioreactor to allow the biology team to conduct experiments whose aim can be, for instance, characterisation of an inducible promoter (pLac in PARSE's case). It is a critical step to determine the levels of growth-inhibiting effect caused by the growth slowing genes, and in the future, the system will gather data in a more efficient and versatile way via its automation.

The hardware modelling was split into different subsystems and then a model was developed joining the subsystems together to explain the integrated reactor system. The subsystems are interconnected between themselves as shown in Figure 0.

Figure 0. Schematic of the hardware model in Simulink, the subsystems are in blocks and their naming is in blue. The connections and the most important interactions and variables produced by the subsystems are represented by the lines. If the line is on the right side of the block, it is an output, and if it is on the left, it is an input to the subsystem. The names of the variables are written in orange.

Theory and Methodology

Growth

The growth of the biomass was treated similarly to what was described in the biological modelling. The growth rate was calculated by:$$\begin{equation}\mu = \frac{\mu_{max}*s}{K_{s}+s}\end{equation}$$ as \(s\) is the mass of substrate present in the vessel, \(K_{s}\) is the half-maximum constant of the process, \(\mu(s)\) is the growth rate and \(\mu_{max}\) is the maximum growth rate reached by the system. The rate of change of the mass of substrate \(\frac{ds}{dt}\) was calculated during the simulation with the equation: $$\begin{equation}\frac{ds}{dt} = -y^{-1}*\mu(s)*x(t) + C_{s0}*v_{in} - \frac{s(t)}{V}*v_{out}\end{equation}$$ this expression is equivalent to an equation shown in the biological modelling. The dilution rate, \(D\), was written as the inflow and outflow of fluid, \(v_{in/out}\) and the concentration of the substrate was changed to \(\frac{s}{V}\), as \(V\) is the volume of fluid in the vessel. The first term of the equation, \(y^{-1}\mu(s)x(t)\), corresponds to the mass of substrate consumed by the biomass. \(x(t)\) is the mass of biomass over time, and \(y\) is the growth yield of the process.

Similarly, the rate of change of the mass of biomass, \(\frac{dx}{dt}\), was calculated in a similar way to another equation in the biological modelling: $$\begin{equation}\frac{dx}{dt} = \mu(s)*x(t) -\frac{x(t)}{V}*v_{out}\end{equation}$$

OD Sensor

The OD sensor was modelled in a simple manner. It considers the concentration of the biomass and whenever this concentration becomes above a certain value, it emits a signal to the fluid control subsystem. The signal initiates the feeding and wash-out process. A small adjustment made to consider a real OD sensor was that it will only detect the concentration if it is above a certain value.

Fluid Control Subsystem

The fluid control subsystem regulates the addition and removal of liquids and provides the total volume of liquid over time in the subsystem. The simulation is adjusted to work as a turbidostat, and the volume inside the vessel varies as shown in Figure 1. When the simulation is run, the control of the subsystem checks to see if the volume in the vessel is equal to the desired volume for the experiment. If the volume is not the same as the desired volume, it adds fluid until it reaches the desired value. Then, no volume is changed until the OD sensor is triggered. Whenever it sends the signal, the control of the subsystem turns the in and outflow to replace a certain percentage of liquid. The in- and outflow are assumed to be equal and constant.

Figure 1. How the volume changes over time in the simulation. In the start of the simulation, it can be seen an increase to the volume for it to achieve the desired value (in this case 15mL). Then, the volume remains unaltered as it is assumed that when the OD sensor is triggered, the inflow is the same as the outflow.

Whenever the control system activates the pump, there will be an in- and outflow. Which means: \(\frac{dV}{dt} = inflow – outflow \). However, as inflow = outflow, the rate of change of volume over time will be 0. The control system is kept activated for a certain period,\(T\), until a percentage of the fluid is replaced. This percentage is set by the user and the period that the pump is activated is calculated by:$$\begin{equation}T = \frac{V(t)*dilution}{inflow}\end{equation}$$ The volume of liquid added in the period \(T\) corresponds to the volume necessary to replace the ratio of liquid aimed for.

The simulation can run infinitely if it is assumed that an infinite amount of fluid is available. However, as there is usually a limited amount of volume that can be added, the user can set a maximum amount of volume for the inflow. After this limit is reached, the control system stops with any in- and outflow. This can be seen in Figure 2 where the maximum amount of volume added is 22.5 mL. In Figure 2, the arrows illustrate a few instances where the OD sensor was activated. In these instances, the same amount of volume is added and removed from the system, which causes the step in both the volume added and removed. At a time of around 450 s, all the lines become constant, with no change. That happens since all the available volume for the inflow was pre-set and the simulation does not change the volume independently.

Figure 2. Graph demonstrating how volume changed over time. The blue line corresponds to the total volume that has been added to the subsystem, and the yellow line corresponds to the total volume removed. The orange line [\(V(t)\)] represents the volume inside the vessel at a certain time.

The fluid control subsystem interacts with the other subsystems. It provides the volume of liquid present in the vessel at a given time, \(V(t)\). And, in the simulation model it indicates to other subsystems when there is fluid added or removed. The latter has the most impact across the thermal and growth subsystems. The effect it has can be seen within each respective subsystem.

With regard to the growth subsystem, it assumes that all cells and substrates are equally distributed in the fluid. Hence, when some liquid is washed out, a certain number of cells and substrate are removed from the solution. Furthermore, when liquid is added to the vessel, substrate is also added to the fluid. The amount of substrate added, \( s_{added} = C_{s0} * inflow\). As \(C_{s0}\) is the concentration of substrate in the solution being added and the inflow is the rate at which fluid is fed into the vessel. In Figure 2. it was shown that there may be a limit of liquid that once fed in would stop the in- and outflow. Since there is no outflow of liquid, there is no biomass removed from the system, which causes the quantity of biomass to increase exponentially over time. Furthermore, as fluid is not added to the subsystem, no more substrate is added either, which causes the biomass to stop growing after the remaining substrate has been consumed. This can be seen in Figure 3. The lines initially are the same, until around 480 s when the maximum volume of inflow is achieved. At this moment, the OD sensor is turned off and the biomass is allowed to grow exponentially again. Opposingly, the infinite supply of fluid continues in the repeating pattern as the OD continues to trigger and there continues to be in and outflow of fluids.

Figure 3. Comparison of the amount of biomass over time when there is no liquid limit (blue line) and when there is a maximum amount of liquid to be added (orange line). After the set volume of liquid is used up, the line that corresponds to the limited inflow of liquid grows exponentially, until it reaches a peak when it consumes all the substrate in the media.

Thermal Subsystem

The thermal subsystem adjusts the temperature of the bioreactor by turning on and off the heating element, and calculating any heat loss from the vessel. The simulation considers the system to be split into the vessel and the fluid, assuming that the temperature of the fluid and the vessel are equally distributed, and any changes to it will happen instantaneously. The vessel is assumed to take a cylindrical shape and there is no heat transfer from the base of the cylinder. The fluid is assumed to have the same physical characteristics as distilled water and the vessel has no irregularities or impurities.

The thermal system is regulated by a control system. Depending on the temperature, it will provide different values of power to heat the system. The thermal subsystem regulates how the power provided to the system in a cycle interacts with the different components. All the power provided must be either absorbed by the liquid, absorbed by the vessel, or lost to the environment. Hence: $$\begin{equation}\sum Power = P_{provided} = P_{liquid} + P_{vessel} + P_{lost}\end{equation}$$

The power provided, \(P_{provided}\), is regulated by a feedback loop within the thermal control system. The heat transferred to the liquid increases the temperature of the fluid following the equation: $$\begin{equation}Q = mc \Delta T\end{equation}$$ as \(Q\) is the heat provided to the fluid, \(m\) is the mass of the fluid and \(\Delta T\)is the change in temperature (final – initial), and \(c\) is the specific heat capacity of the fluid. Since the energy provided \(Q(t)\) is dependent on time, and the final temperature, \(T(t)\), will also depend on time, the equation can be rearranged to: $$\begin{equation}\frac{dT}{dt} = \frac{1}{mc}*\frac{dQ}{dt} = \frac{P}{mc}\end{equation}$$ \(\frac{dQ}{dt}\) is the rate of change of energy, or the power, \(P\), provided to the system. However, in the model, all power added into the system is transferred to the vessel. The vessel exchanges temperature with the fluid and the power provided to the fluid has value: $$\begin{equation}P_{liquid} = h* SA (T_{v}-T_{l})\end{equation}$$ as \(h\) is the thermal conductivity of water, \(SA\) is the surface area of contact between the vessel and the liquid and \(T_{v} – T_{l}\) is the difference in temperature between the vessel and the liquid. The surface area in contact with fluid has been assumed to be solely the circular base of the cylinder. Hence, \(SA = 2*\pi*r*height\), as \(r\) is the radius of the cylinder. The value of \(h\) varies with the temperature of the water. However, the change in the value is not meaningful; it varies from \((1.095 - 1.142) *10^{-4}\)an increase of 4%. Hence, this model assumes that this value will remain constant. To understand how to calculate the value of \(h\), see [1]. This way, the value for \(P_{liquid}\) is calculated.

If the round sides of the cylinder are thought to be slabs, the temperature of the inside of the cylinder, which touches the liquid, should be the same temperature as the liquid. Similarly, the outer side of the cylinder should have a similar temperature to the environment it is in. There will be a change in temperature over the width of the cylinder illustrated in Figure 4. an illustration of how the temperatures should be distributed across the wall. Equation (9) can be used to equate the power transferred, \(P_{lost}\), from the fluid to the environment [2]: $$\begin{equation}P_{lost} = \frac{-k*2\pi r*height}{d}*(T_{l} - T_{environment})\end{equation}$$

Figure 4. Sketch of how the temperatures of the walls of the figure should be distributed. The areas in blue should have a temperature similar to the fluid. The areas in red should have a temperature that approaches the external temperature the closer it is to the outside.

In (9), \(k\) is the thermal conductivity of the material that the cylinder is made of, \(r\) and \(h\) are, respectively, the radius and height of the cylinder and \(d\) is the width of the walls of the cylinder.

After \(P_{liquid}\) and \(P_{lost}\) are calculated by the simulation and \(P_{provided}\) is given by the control system, (5) is used to calculate the power absorbed by the vessel. After the power is calculated, (4) is used to calculate the change in temperature of the vessel and the change in temperature of the fluid. This rate of change is passed through an integrator block in Simulink and the temperature \(T(t)\) of the fluid and the vessel are changed accordingly.

The model allows the user to choose the option to have insulation. The addition of insulation does not change the main mechanics of the subsystem. However, a few adjustments are required. For instance, the power lost to the environment needs to account for the decreased heat transfer through insulation. For this case, insulation is assumed to be connected to the material of the vessel (illustrated in Figure 5.), with a width of \(d_{ins}\). To account for the thermal conductivity of both the tube and the insulation, the equation is corrected: $$\begin{equation}P_{loss} = \frac{T_{l} - T_{environment}}{SA*(\frac{d_{ins}}{k_{ins}}+\frac{d_{vessel}}{k_{vessel}})}\end{equation}$$ as \(k\) is the thermal conductivity of the material and \(d\) is the width of the material. Both \(k\) and \(d\) can be of the insulation or the vessel depending on the subscript, ins, or vessel, respectively. \(SA\) is the surface area where the power loss happens, and \(T\) is the temperature, either of the fluid or the environment, depending on the subscript. The insulation is assumed to not receive any heat from the power of the vessel as it is supposed to have a small heat transfer. Hence, it is not considered when calculating the increase in temperature of the reactor vessel.

Figure 5. Illustration of the connection between the insulation and the test tube (vessel). The insulation is represented by the orange region and the test tube/vessel is represented by the black region.

The addition of insulation influences the temperature as confirmed by Figure 6 (a) and (b). In (a), there is a comparison of the temperature of the liquid with and without insulation. Both the simulations were run with the same parameters and intended to reach a temperature of 37 °C. However, as the blue line (without insulation) indicates, the desired temperature was not reached, and it stabilised at around 30 °C. The fact that the temperature in the simulation without insulation could not reach the desired temperature is explained in (b). In (b), the power loss over time was compared between the simulations. In the simulation without insulation, the blue line, the power lost exponentially increased until it reached the maximum wattage being provided to the system. This caused all the heat added to the system to be lost and the temperature of the fluid would not increase. Hence, the simulation demonstrates that for desired temperatures, the presence of insulation is essential for the experiments.

Figure 6. (a), (b). The figures compare simulations run with and without insulation. On the left, the temperature is compared while on the right, the power loss is compared. The blue line corresponds to the simulation without insulation and the orange line corresponds to the simulation with insulation.

The thermal system is dependent on the fluid control system in the simulation. Although the volume of fluid does not change over time, the new fluid added to the system is assumed to have a lower temperature than the fluid inside of the system. To calculate the equilibrium temperature between the old and new fluids, the following equation is used: $$\begin{equation}Q = m_{old}*c*(T_{f} - T_{old}) +m_{new}*c*(T_{f}-T_{new})\end{equation}$$ as \(m\) is the mass of fluid, \(c\) is the specific heat density. \(T_{f}\) is the average temperature between the old and new media and \(T_{old/new}\) corresponds to the temperature before the new media was introduced. In the simulation, \(T_{new}\) is the same as the temperature of the environment and \(T_{old}\) corresponds to the temperature of the fluid after it goes through the integrator. Since this interaction presents only heat transfer between the fluids, \(Q = 0\). If the specific heat capacity, \(c\), is assumed to be the same for both fluids, \(c\) can be canceled, which allows (11) to be rearranged: $$\begin{equation}T_{liquid} = \frac{T_{old}*m_{old}+T_{new}*m_{new}}{m{old}+m_{new}}\end{equation}$$ \(T_{liquid}\) is then considered to be the temperature of the fluid for the next iteration of the process.

Agitation Subsystem

The agitation subsystem only performs a single function: activate a motor to spin a magnetic stirrer. The control subsystem turns on whenever the simulation starts, and turns off whenever the simulation ends. However, although the agitation subsystem can be seen as “symbolic” at this stage, it calculates the mechanical power transferred to the liquid during the simulation. The subsystem calculates the equation: $$\begin{equation}P = k_{n} *\rho * \omega^{3} *d_{imp}^{5}\end{equation}$$ As \(P\) represents the power transferred to the liquid, \(k_{n}\) the power number (Newton number) of the stirrer, \(\rho\) the density of the fluid, \(\omega\) the angular velocity of the stirrer, and \(d_{imp}\) the diameter of the stirrer. The density of the fluid can be adjusted in the simulation, however, all simulations performed for this report use \(\rho = \rho_{water}\) for density.

Aeration Subsystem

In small-scale experiments, the amount of dissolved oxygen (\(dO_{2})\) is commonly assumed to be sufficient for any culture growth. However, for large-scale experiments and bioreactors, this assumption most times is not true. Hence, the subsystem was developed to complete the system model and to demonstrate that in small-scale experiments the amount of oxygen at the start is sufficient. The subsystem uses the assumption that there is a constant inflow of sterile air into the fluid and the whole surface area can dissolve oxygen. Furthermore, the gas inflow is assumed to be thermally isolated from the liquid, which means that there is no heat lost or transferred to the gases.

The rate of change on the concentration of \(dO_{2}\) being dissolved is given by: $$\begin{equation}\frac{dC_{dO_{2}}}{dt} = k_{l}a*(C_{sat}-C(t))\end{equation}$$ as \(C_{sat}\) is the concentration of saturated oxygen in a certain temperature, \(k_{l}a\) is the volumetric oxygen mass transfer coefficient of the experiment and \(C(t)\) is the actual concentration of \(dO_{2}\) at a time step. The \(k_{l}a\) is given by (15) where \(C\), \(\alpha\), \(\beta\), and \(\delta\) are constants, \(v_{s}\) is the superficial velocity of the gas, \(V\) is the volume of the liquid, \(\mu\) is the viscosity of the liquid and \(P_{g}\) is the power of the gas. $$\begin{equation}k_{l}a = C*v_{s}^{\alpha}*(\frac{P_{g}}{V})^{\beta}*\mu^{\delta}\end{equation}$$

The saturated concentration of \(dO_{2}\) varies with temperature. The relationship between these quantities is complex. However, for the range of temperature considered in the model (15 – 40 °C), it assumes that the saturation concentration of \(dO_{2}\) varies linearly with temperature. Figure 7. represents the real relationships of temperature and the concentration in blue and the simplified linear relationship of the variables in orange. The relationship found performing a linear regression with the data was: $$\begin{equation}C_{sat} = T*(-0.145)+12.0\end{equation}$$ the linear regression provided a value for the correlation coefficient of \(R^{2}= 0.989\). As the \(R^{2}\) value of the relationship is close to 1, it can be said that the relationship is strong. This indicates further that this simplified linear relationship can be used.

In (15) the variables \(P_{g}\) and \(v_{s}\) are calculated as referenced in [3]. The value of the superficial velocity, \(v_{s}\), was calculated by \(v_{s} = \frac{airflow}{cross-sectional area}\). Since the bioreactor is modelled as a cylinder, the cross-sectional area is going to be a circle with area \(A = \pi r^2\). The value of \(P_{g}\) is calculated based on the equation: $$\begin{equation}P_{g} = \frac{P}{\sqrt{1+\frac{Y*v_{s}}{\sqrt{g*d_{imp}}}}}\end{equation}$$ as \(P\) is the mechanical power calculated in the agitation subsystem, \(Y\) is a constant, \(g\) is the acceleration due to gravity and \(d_{imp}\) the diameter of the stirrer. According to [3], the constants in equation (15) have the value \(C =0.026\), \(\alpha = 0.4\), \(\beta = 0.5\), \(\delta = 0.55\). In (17), \(Y = 750\) for one stirrer.

Figure 7. Comparison between two models of how the solubility of oxygen changes with temperature. The blue line indicates the ”real relationship” between the variables, with all the data points collected being considered. The orange line corresponds to the simplified version where a linear regression was performed to indicate the correlation between the variables.

Equation (17) considers the oxygen released/attached to the fluid inside the bioreactor. Another aspect that was considered in the model, was the oxygen consumption of the biomass. A simple proportional relationship was considered where the amount of oxygen consumed is equal to the mass of biomass multiplied by a constant. This constant was the mass of oxygen consumed by one gram of biomass. Hence, with the information of the mass of oxygen consumed by the biomass, the final equation was reached: $$\begin{equation}\frac{dC_{dO_{2}}}{dt} = k_{l}a * (C_{sat}-C(t)) - x* B\end{equation}$$ as \(x\) is the mass of biomass present in the solution and \(B\) is the constant, mass of oxygen consumed by a gram of biomass.

Results and Analysis

System Modelling and Simulation

The Simulink model and the supporting documentation can be found in the iGEM Sheffield GitLab repository.

The controlled variables calculated by the Simulink model were gathered into two different scopes. They were separated into the variables that were considered the “Physical” and “Biological” parts of the simulation. These are presented in Figures 8., 9., respectively.

Figure 8. Demonstrates the physical aspects of the simulation. The volume of fluid over time, the concentration of \(dO_{2}\) during the simulation, and the temperature of the vessel and the liquid over the simulation are shown, respectively. In the volume vs time graph, after the initial injection of liquid, the volume remains constant. Although the value is shown as a constant, as it was mentioned before, there is an inflow/outflow of volume happening at different points. In the \(dO_{2}\) concentration, the exponential decrease is due to the temperature dependency of this value. It stabilises after some time when the temperature stabilises. The simulation in this case was run with sufficient insulation that the power loss would not be greater than the power provided. This has caused the temperature of both the vessel and the fluid to stabilise after it reaches the expected temperature of 37 °C.

Figure 8. Graphical representation of the Physical variables produced by the Simulink model. On the left, the volume of fluid over time is presented. In the middle, the concentration of dO2 during the simulation is shown. On the right, the temperature of the vessel is demonstrated by the yellow line and the temperature of the liquid is seen in the red line.

The graphs presented in Figure 9. correspond to the biological aspect calculated by the simulation. The left graph shows how the concentration of biomass (\(Cx)\) varies over time. The peaks and droughts in the middle of the graph correspond to the feeding/wash-out process after the OD sensor is triggered. The exponential increase indicates that all the volume has been used and there is no extra volume to be added. The two following graphs present a similar trend. Whenever there is an increase in the mass of the substrate, there is an increase in the growth rate. These increases correspond to the addition of substrate-rich solution to the vessel, which increases the mass of the substrate.

Figure 9. Graphical representation of the biological variables produced by the Simulink model. On the left, the concentration of biomass over time is presented. In the middle, the mass of the substrate during the simulation is shown. On the right, the growth rate of the cells during the simulation is shown.

System Design and Engineering

iGEM Sheffield's System Design and Engineering adheres to a systematic and structured approach for conceptualising, planning, and developing a bioreactor that aligns with defined objectives, functionality, and performance criteria. This comprehensive effort involves designing, integrating, and optimising diverse subsystems and their components to create an efficient, reliable, and cost-effective solution [4].

The analytical results yielded from the Simulink-based modelling and simulation of a turbidostat reactor system have been used throughout the System Design and Engineering process to inform and analyse decisions regarding the integration of mechanical and electronic components, subsystems, and the programming of control algorithms. For instance, the aeration subsystem in the model demonstrated that it was not necessary to be preoccupied with the amount of \(dO_2\).

Objective

iGEM Sheffield aimed to develop a modular, multiplexed turbidostat explicitly built for growth characterisation under dynamic conditions. This system is intended for growth modulation experiments across various scales, offering a more streamlined and cost-effective alternative to intricate and expensive microplate spectrophotometers. This type of bioreactor makes it possible to automatically build an inducer vs OD standard curve, and to obtain several similar curves for growth characterisation purposes in parallel. With its several cells that may contain different microbial cultures, multiple cultivation experiments under varying conditions can be run simultaneously and analysed in real-time while maintaining turbidity at a constant level.

Growth characterisation is the process of studying and characterising how microbial populations grow under specific conditions through controlled cultivation. By measuring parameters such as growth rate, lag phase, stationary phase, biomass yield, and other relevant factors, researchers can gain insights into the growth behaviour of microorganisms [5].

Engineering Design Process

The development of the bioreactor system followed an interdisciplinary, iterative engineering design process characterised by a structured yet adaptive approach. This method involves constructing physical and virtual models in incremental stages, each representing varying levels of fidelity, facilitating testing, feedback, and optimisation [6]. The iterative engineering design approach helped mitigate the impact of changing and evolving requirements.

Figure 10. Interdisciplinary Hardware System Modelling and Engineering team. From left to right: Berardo Manuel Sanchez Tafolla - PhD Student (3° Year), PhD Biomedical Science Development, Regeneration and Neurophysiology, School of Biosciences; Lucca Johann Leal - BSc Student (2° Year) Theoretical Physics, Department of Physics and Astronomy, Faculty of Science; William Nathaneal Santoso - BEng Student (2° Year) Mechatronic and Robotic Engineering, Department of Automatic Control and Systems Engineering, Faculty of Engineering. University of Sheffield — 2023

Design Research

In the build-up to the system design and engineering phase, required and preferred specifications for a small-scale multiplexed bioreactor were explored — including considerations about usability, laboratory applicability, materials and technologies, sustainability and environmental impact. Further insights into the design of bioreactors were attained through interviews conducted with academic experts in the Faculty of Engineering at the University of Sheffield, particularly concerning fluid control systems and the development of a multispectral fluorometer.

Figure 11. Educational Bioreactors in the Bioengineering Laboratory, Diamond Engineering Building, Faculty of Engineering. University of Sheffield — 2023

Design Iteration I

The first design iteration focused on the development of a modular, small-scale bioreactor for cultivating and harvesting multiple microbial biomasses at a reactor cell volume of 150 - 200 mL. This initial bioreactor design was envisioned to be easily customisable and to have its individual cells operate in parallel or series as either turbidostat, chemostat or fluorostat. For preliminary review, the reactor was prototyped as a low-fidelity (LoFi) physical model—using cardboard and synthetic polymers—to gain insight into the geometrical form and approximated dimensions, assumed logic for the control engineering, potential component fitting, and assembly complexity.

Figure 12. Representation of the first design iteration for a small-scale bioreactor for cultivating and harvesting multiple microbial biomasses at a reactor cell volume of 150 - 200 mL. Left: the multiplexed bioreactor as a 2D top-down concept graphic. Right-top: A single reactor cell as cardboard and polymer prototype. Right-bottom: Two reactor cells as cardboard and polymer prototypes to exemplify module connection.

This initial design concept was discarded to accommodate changes in the direction and focus of the biological experiments conducted by iGEM Sheffield 2023 within the context of growth modulation.

Design Iteration II

For the second design iteration, learnings from the previous prototype were extracted and used as a baseline for the general functionality of a multiplexed bioreactor, particularly on reactor control engineering and sensor design and deployment.

In addition, further analysis of the adjusted biological processing applied in iGEM Sheffield's experimental design was conducted to identify areas for improvement. As a result the workflow of growth characterisation—usually requiring the use of microplate spectrophotometers—became the targeted use case for the second design iteration. Hence, the new reactor system was conceptualised to operate as a multiplexed turbidostat for cultivation, real-time analysis, and study of multiple microbial cultures at a sample-sized reactor cell volume of maximum 25 mL.

To expedite the review process, physical models and representations were prototyped as low-fidelity 2D models electronic control system schematics and breadboard assemblies. The review process for the second design iteration included interviews with PhD students of the School of Biosciences at the University of Sheffield, particularly with regard to requirements and specifications for growth characterisation in biological experiments and laboratory settings.

Figure 13. Review of the requirements and specifications for growth characterisation in biological experiments and laboratory settings with PhD students of the School of Biosciences. University Sheffield — 2023
Figure 14. Left: Draft concept schematic outlining the reactor cell’s design structure and component integration. Right: Draft concept schematic outlining the configurable modular design of the multiplexed turbidostat
Figure 15. Draft concept schematics outlining the electronics layout and power supply distribution for the multiplexed turbidostat.
Figure 16. Breadboard assembly of the basic circuit design for the reactor cell components in the Electronic and Electrical Engineering Laboratory, Diamond Engineering Building, Faculty of Engineering. University of Sheffield — 2023

Design Iteration III

Building on the validated concepts of the second design iteration, the third iteration focused on developing mechanical and electronic specifications for production. Conducting thorough examination and design from the sub-system level to the integrated system ensures that each component, from mechanical elements to electronic circuitry, is planned and designed to align precisely with the intended objectives, design specifications and operational requirements

Electronic and Electrical System Design

The electronics of the multiplexed turbidostat system are segmented into two distinct circuit types, interconnected through magnetic connectors: the control unit circuit and the reactor cell circuit. The control unit circuit is primarily responsible for the power management and the fluid control subsystem, whereas the reactor cell circuit mainly oversees the thermal subsystem, agitation subsystem, and sensor readings.

Power supply for the multiplexed turbidostat relies solely on a wall adapter, which steps down mains voltage to 12V within the control unit circuit. Subsequently, this 12V power supply is distributed to the adjacent reactor cell circuits. This centralised power distribution approach affords the system greater flexibility in the arrangement of cells, without the need for multiple bulky power supply units.

The following sections provide a comprehensive overview of the specific components housed within each type of circuit, along with their respective functionalities within different subsystems.

Control Unit
Figure 17. Full electronics schematic for the fluid control unit of multiplexed turbidostat. Note: The fluid control unit schematic only depicts connections to one set of stepper and servo motors, not the planned three, in order to avoid repetition.

Power and Control System

Fluid Control Subsystem

Reactor Cell
Figure 18. Full electronics schematic for the reactor cell of the multiplexed turbidostat.

Power and Control System

Agitation Subsystem
Thermal Subsystem

However, it's important to note that the PTC heating element, by default, only self-regulates to its rated temperature when powered. To precisely control and maintain the desired temperature, an external relay or MOSFET is required to execute Proportional-Integral-Derivative (PID) control. In this capacity, the chosen component for temperature control is the Gravity MOSFET Power Controller. Its selection was based on its compatibility with the voltage and current requirements of the PTC heating element, rendering it a suitable choice for achieving the precise temperature control needed for the system.

Embedded Software Architecture

The embedded software consists of a main source file and series of helper files to carry out specific functions. The software will be written in MicroPython. Each subsystem will have their own corresponding helper files to define their configurations and functions that will be executed in the main source file. External libraries for each component are also used to simplify communication and interfacing with said components. This is done to facilitate abstraction in the main source files, allowing them to be cleaner and more readable. Breaking down the software in this way also helps with debugging and troubleshooting, as it allows the root cause of any given error to be more easily traced.

3D Printing Materials

The holding structures, sleeves, and inlays will be fabricated through 3D printing utilising two different filament materials: Polyethylene Terephthalate Glycol (PETG) and Polylactic Acid (PLA). PETG will be the material of choice for components situated in close proximity to the heating element. This choice is due to PETG's high mechanical strength and low thermal conductivity, mitigating the risk of melting or warping due to the heat exerted by the PTC heating element. PLA, by contrast, will be employed for the remaining structures primarily due to its cost-effectiveness. The strength characteristics of PETG are not necessary for these components, rendering PLA the better option for this purpose.

Table of Components

Component

Quantity

Power Supply

Power Supply Adapter

1

5V Buck Converter

7

Magnetic Power Connectors (5 Pin)

6

Reactor Cell (4)

Raspberry Pi Pico W

4

Noctua NF-A4x10 PWM

4

Neodymium Magnets 

8

PTFE Stir Bar

4

PTC Heating Element

4

MOSFET Power Controller

4

Holding Structure (incl. Inlays)

4

Light Emitting Diodes (LED) - Power & State

8

Gravity MLX90614-DCC IR Temperature Sensor

4

Light Emitting Diode (LED) - OD

4

Light Emitting Diode (LED) - Fluorescence

4

Phototransistor - OD

4

Phototransistor - Fluorescence

4

Resistors 

20

Emission Filter

4

Excitation Filter

4

Glass Specimen Vials (20 mL)

4

Vial Screw Cap

4

Vial Screw Cap Sleeve

4

Steel Capillary Tubes

12

Vial Holding Sleeve

4

Control Unit

Raspberry Pi Pico W

1

Glass Media Vials (100 mL)

2

Glass Media Vials Screw Caps

2

MG996R 360° (Metal Gear) Servo Motor

3

MXT Pinch Valve

3

Peristaltic Pump Head

3

NEMA 17 Stepper Motor

3

TB6612 Motor Driver

3

Silicon Tubing

12

Holding Structure (incl. Inlays)

1

Resistors

2

Light Emitting Diodes (LED) - Power & State

2

General Items

HATCHBOX PETG Filament - 1 kg

1

PLA Filament - 1 kg

1

18 AWG Wiring - 5.5 m

1

Jumper Wires - 120 pcs

1

Table of components for the multiplexed turbidostat assembly — featuring one fluid control unit and four reactor cells.Please Note: This list offers an overview of key components based on the outcomes of Design Iterations I to III. The final component list may undergo revisions following further adjustments due to testing of high-fidelity prototypes for critical subsystems prior to production.

The third design iteration is concluded with high-fidelity prototyping and testing of critical subsystems, such as the custom design for a multiplexed pinch valve.

The multiplexed pinch valve is employed in the turbidostat reactor as part of the fluid control system to regulate the flow of culture media, nutrients, or other agents between the control unit and the reactor cells. By controlling the pinch valves, researchers can adjust the flow rates (in- and outflow) to maintain specific culture conditions, such as constant optical density or cell growth rate, simultaneously in multiple cultures or experiments.

The goal of the high-fidelity prototyping is the design of detailed, functionally accurate representations of critical reactor system components before final production. They include both 3D CAD models and physical builds. These high-fidelity prototypes closely mimic the turbidostat’s actual design, component integration, and functionality, providing a platform for rigorous testing and validation. This approach helps identify and address potential issues and refine system performance while ensuring the multiplexed turbidostat meets its intended design specifications and operational requirements.

Example: High Fidelity Prototyping of Multiplexed Pinch Valve
Figure 19. Static 3D model of the multiplexed pinch valve, employed in the turbidostat bioreactor to regulate the in- and outflow of fluids for up to four reactor cell vessels.
Video 1. 3D exploded view animation of the multiplexed pinch valve, employed in the turbidostat bioreactor to regulate the in- and outflow of fluids for up to four reactor cell vessels.
Figure 20. Images of the 3D printed, fully functional high fidelity prototype of the multiplexed pinch valve (Fordham-Brown concept).
Multiplexed Pinch Valve Pressure Integrity and Leakage Test

To confirm a secure seal, the closed multiplexed pinch valve underwent a pressure integrity and leakage test with fluid pressurised to 40 psi, while submerging the pipe's end in water to detect any escaping air bubbles. The pressure gauge recorded the test conditions.

Video 2. Pressure integrity and leakage test conducted to confirm a secure seal of the multiplexed pinch valve when closed.

Conclusion

Building an open source, multiplexed turbidostat system requires to a greater degree original, purpose-built engineering concepts for critical parts such as the multiplexed pinch valve or the integration of the reactor cell sensors. This is necessary to simplify the sourcing of required complex or otherwise costly components and improve maintainability and operation while not compromising the system’s reliability and efficiency. These original engineering concepts increase the time-to-build due to extensive research-driven prototyping stages to identify and test different materials, components, and design structures.

To accommodate research-driven prototyping, iterative design in the engineering of the modular, multiplexed turbidostat has proven to be an effective approach to account for changing or evolving requirements. It has also allowed for the consultation of experts from different disciplines at different stages of the prototyping without impeding the design process.

Further work towards the manufacturing stage of the engineering process will be informed by the design and testing of high-fidelity prototypes for critical reactor system components such as the modular, multiplexed pinch valve.

Value Proposition

The modular multiplexed turbidostat is being designed as an alternative to conventional microplate spectrophotometers in specific biological research applications by offering several operational advantages and cost-saving benefits:

However, microplate spectrophotometers may still be preferable for assays that require greater precision in optical measurements, specific assay compatibility, high sample throughput, or specialised filters and optics not provided by a multiplexed turbidostat.

To further enhance the modular, multiplexed turbidostat, it has been designed with several IoT (Internet-of-Things) redundancies and capabilities, providing following additional advantages in its deployment:

Introducing robust IoT redundancies and capabilities into the modular multiplex turbidostat enhances its deployment. With features such as remote monitoring and control, automated data collection and storage, and seamless integration into a wireless turbidostat network, the system provides researchers and educators with greater flexibility, efficiency and reliability.

Future Work

System Modelling and Simulation

For future work on the modelling and simulation of the bioreactor the complexity of the agitation and aeration subsystems could be increased. In the current model, these systems are underdeveloped as there was no large incentive to increase their complexity due to the design specifications for the 2023 iGEM Sheffield hardware system. However, for larger-scale bioreactors these systems are more complex as the production of carbon dioxide is higher. The greater increase in dissolved carbon dioxide could cause issues with the growth of cells [8]. Hence, the addition of a control system to the agitation and aeration subsystem to regulate the airflow and speed of the rotor would be an essential to improve this model for the design of bioreactors with a larger working volume.

Another aspect that could be further developed is the integration of the different steps of the modelling performed. The integrated model could work as follows: the data provided by a SimBiology model would be transferred to a reaction process script, which would analyse it and provide the necessary parameters into the Simulink model. The results provided from the Simulink model could then be used in a Simscape model that would be a representation of the physical system used by the system's design team to make design decisions.

System Design and Engineering

Multispectral Fluorometer Design

The attempted design of a compact, low-cost multispectral fluorometer with laboratory standard accuracy has posed several engineering challenges. Further research and analysis is required to evaluate component choices and designs for beam splitting, excitation and emission filtering, and a solid-state or photodiode detector array to determine the feasiblility of such a concept.

Manufacturing

Further testing of the sub-systems—in particular its electronic design and control programming—is recommended before transitioning the modular, multiplexed turbidostat to the manufacturing stage. The range of tests will include functional testing, performance testing, integration testing, electrical safety testing, signal integrity testing, and Failure Mode and Effects Analysis. The availability of high-fidelity prototypes of the critical subsystems, and their test results, will provide the groundwork for manufacturing the multiplexed turbidostat according to the required design specifications and operational requirements while reducing the likelihood of redesigns or modifications during production.

To field-test the multiplexed turbidostat in a pilot phase, it is planned to manufacture a fully operational system — including a control unit and up to four reactor cells. The manufacturing will take place in the iForge Makerspace and Electronic and Electrical Engineering facilities at the University of Sheffield. It will involve 3D printing, electronic and mechanical component assembly, embedded control programming, the integration of the subsystems, quality testing and design reviews. Further learnings derived from the pilot phase can then be used to inform future iGEM Sheffield bioreactor designs.

Hardware References

[1] M. Bahrami, “Forced Convection Heat Transfer,” 2011. Accessed: Sep. 30, 2023. [Online]. Available: https://www.sfu.ca/~mbahrami/ENSC%20388/Notes/Forced%20Convection.pdf
[2] B. Reisfeld, “Heat Transfer — Introduction to Chemical and Biological Engineering,” www.engr.colostate.edu, 2023. https://www.engr.colostate.edu/CBE101/topics/heat_transfer.html (accessed Sep. 30, 2023).
[3] J. Vanags and A. Suleiko, “What is k L a in bioreactors?,” Bioreactors.net, 2021. Accessed: Sep. 14, 2023. [Online]. Available: https://741b25f1-5d67-443e-98a9-d8a2a8564dab.filesusr.com/ugd/7abc3e_10bd1c7f1a494ad88f009673bfebdfd1.pdf?index=true
[4] D. M. Buede and W. D. Miller, The engineering design of systems : models and methods. Hoboken, New Jersey: Wiley, 2016.
[5] R. M. Maier, C. P. Gerba, and I. L. Pepper, Environmental microbiology, Second edition. San Diego: Academic Press, 2009, pp. 37–54.
[6] R. Adams, Understanding design iteration: representations from an empirical study, in Durling, D.and Shackleton, J. (eds.), Common Ground - DRS International Conference 2002 , 5-7 September, London, United Kingdom, 2002
[7] P. Mira, P. Yeh, and B. G. Hall, Estimating microbial population data from optical density, PLOS ONE, vol. 17, no. 10, p. e0276040, Oct. 2022, doi: https://doi.org/10.1371/journal.pone.0276040.
[8] D. R. Gray, S. Chen, W. Howarth, D. Inlow, and B. L. Maiorella, “CO2 in large-scale and high-density CHO cell perfusion culture,” Cytotechnology, vol. 22, no. 1–3, pp. 65–78, 1996, doi: https://doi.org/10.1007/bf00353925.

Biological Modelling

Introduction

The building of genetic constructs encoded on bacterial plasmids has represented a fundamental strategy for the study of gene function in the molecular and synthetic biology fields. A typical transcriptional unit in prokaryotes contains a promoter (the binding region for the RNA Polymerase), a ribosome binding site (RBS), a coding sequence (gene) and a terminator (the unbinding site of the transcription complex). Depending on its activation dynamics, a promoter can be classified as constitutive, which is always able to initiate mRNA transcription or inducible, when it requires the presence of an external molecule (typically a small molecule inducer and its associated transcription factor protein) to modulate the formation of the transcription complex.

When testing the function of a heterologous plasmid and its contained genes, two critical factors must be considered to ensure the expected function in the cell: an uninterrupted expression of the transcriptional unit at the right times and amounts; and the absence of undesired effects of this transcription that may otherwise compromise normal cell function. With this in mind, fluorescent reporter genes have become a valuable tool as they enable the analysis of transcriptional unit expression by the production of a single molecule of a fluorescent protein with each round of transcription. Also, their sensitivity, stability and variety of optical properties have earned them reliability and provide them with versatility for these applications.

To test the expression and function of our genetic system and aim to develop a model to describe its expression, we created a computational simulation of how a contained vsfGFP gene would be expressed as a part of our system in optimal conditions. Several parameters that dictate the mechanisms of the involved reactions were considered. Both constitutive and IPTG inducible promoters were included, and their promoter strength metrics were calculated. With this information, the actual transcription level of our developed set of plasmids could be estimated in each condition. On the other hand, the resulting values were compared with the experimental data to test the validity of our initial model.

Method

Determining promoter strength allows the characterisation of inducible promoters, and within PARSE, enables pLac strength to be determined in relation to the Anderson scale. This comparison was completed using RStudio to enable us to generate plots for easier visualisation and calculation of the promoter strength from constructs of pLac- vsfGFP in pET28 and pLac-GS in pET28. Promoter activity (P) is the amount of gene transcription into mRNA. This can also be defined as the polymerases that bind to the promoter and process the entire gene sequence [1], or the amount of polypeptides produced by the process [2]. Using fluorescence to represent the quantity of vsfGFP confines our definition to the amount of vsfGFP. In reality, nascent vsfGFP polypeptides will then undergo a crucial post-translational modification process to produce an active vsfGFP, or - in the case of aberrant expression or folding - undergo degradation by quality control proteases.

However, our model, illustrated in Figure 1 below, makes the assumption that degradation of vsfGFP (D), active or inactive, is 0 (D=0). This means it’s assumed that all vsfGFP produced by expression of the gene will fluoresce. This assumption is made as vsfGFP is closely related to sfGFP; this is a GFP with enhanced stability, and a higher thermostability to standard GFPs [3].

Figure 1: The modelling structure of the protein expression of fluorescent GFP. The rate determining step is assumed to be the translation and transcription of \(GFP_{f}\) which is generalised as P as the promoter activity quantified as \(GFP_{n}\) per cell per hour.

To determine the promoter activity, Equations \(1\) and \(2\) below are first used to describe the general reaction pathway of inactive and active vsfGFP, as illustrated from Figure 1 above. $$\begin{equation}\frac{\partial n}{\partial t}= P- m\cdot n - \mu \cdot n -D_n\end{equation}$$ Where \(P\) is the promoter activity \(GFP_{n}\) per cell per hour, \(m\) is the maturation constant \(h_{-1}\), \(n\) is the amount of inactive GFP (yet to be folded), \(\mu\) is the specific growth rate, and \(D_n\) is the degradation rate of inactive GFP. $$\begin{equation}\frac{\partial f}{\partial t}= m \cdot n - \mu \cdot f -D_f\end{equation}$$ Where \(f\) is the amount of active, fluorescing GFP, and \(D_f\) is the degradation rate of active GFP.

The function of m \(\cdot n\) quantifies the amount of inactive vsfGFP that will be converted into active GFP in a given time.\m represents the maturation constant \(h_{-1}\) which represents the temporal delay between vsf GFP expression on the onset of active fluorescent protein generation in the cells,. \(\mu\cdot n\) or \(\mu\cdot f\) quantifies the amount of inactive or active GFP produced respectively that is dependent on the specific growth rate of the microbial culture. In reference to the Simbiolgy simulation, when the growth rate becomes greater, the accumulation of the viable bacteria is less. This later reduces the amount of inactive and active vsfGFP, making dilution to be among the reaction pathway that minimises the accumulation of fluorescent GFP.

Within a steady state, the amount of inactive vsfGFP can be determined by the equation below [5]: $$\begin{equation}n_{ss}=\frac{P\:-\:D_{n_{\mathrm{ss}}}}{m\:+\:\mu\:}\end{equation}$$ Where \(n_{ss}\) is the amount of inactive vsfGFP at steady state, and \(D_{n_{\mathrm{ss}}}\) is the degradation rate of inactive vsfGFP at steady state.

If we assume that \(D_{n_{\mathrm{ss}}}\) and \(D_{f_{\mathrm{ss}}}=0\) as degradation of vsfGFP is very slow due to its high stability compared to other GFPs.
Then equation 3 can be assumed to be: $$\begin{equation}n_{ss}= \frac{P}{m+ \mu}\end{equation}$$ This can also be determined for the amount of active vsfGFP at steady state. $$\begin{equation}f_{\mathrm{ss}}=\frac{m \cdot n_{\mathrm{ss}}-D_{f_{\mathrm{ss}}}}{\mu}\end{equation}$$ Or $$\begin{equation}f_{\mathrm{ss}}=\frac{m \cdot n_{\mathrm{ss}}}{\mu}\end{equation}$$ for a system that assumes \(D=0\). Where \(f_{ss}\) is the amount of active GFP at steady state, and \(D_{f_{\mathrm{ss}}}\) is the degradation rate of active vsfGFP at steady state.

The relationship of fluorescence for a given time illustrated from F vs t plots from experimental Results is expressed in Equation \(5\) below. $$\begin{equation}\frac{\partial F}{\partial t}= m \cdot N -OD -D_f\end{equation}$$ Where \(F\) is the amount of fluorescence, \(N\) is the amount of non-fluorescence and OD is optical density. As the amount of non-fluorescence increases, the fluorescence intensity also increases however, with the decreasing of the specific growth rate.

The trend of the optical density \(OD\) with time collected from experimental runs is plotted to visualise the growth behaviour of the bacteria. Given the equation: $$\begin{equation}\frac{\partial O D}{\partial t}=\mu \cdot O D\end{equation}$$ Where \(OD\) is the optical density, \mu as the specific growth rate \(h^{-1}\), and t as time \(h\).

The plot is further linearised to \(log_{2} OD\) versus time to determine the specific growth rate \(h^{-1}\) that is derived from Equations (10)-(13).

When both equations are combined, \(\frac{\partial F}{\partial O D}\) can be determined, and this correlates to the amount of fluorescence produced by a single cell as shown in Equation 7 below. $$\begin{equation}\frac{\partial F}{\partial O D} =\frac{\left(m \cdot N-O D \cdot D_f\right) \cdot \partial t}{(\mu \cdot O D) \cdot \partial t} \\=\frac{m \cdot \frac{N}{O D}-D_f}{\mu} \\ =\frac{m \cdot n-D_f}{\mu}\end{equation}$$

Taking that \(f_{ss}\) is equal to \( \frac{\partial F}{\partial O D}\), \(n_{ss}\) is substituted from Equation 3 to Equation 4 where the promoter activity can be determined with known \(f_{ss}\), specific growth rate calculated and maturation constant. Due to the unknown maturation constant of vsfGFP, it is estimated to be ln2 divided by the maturation time of sfGFP which is 13.6 min [6].

$$\begin{equation}P=f_{s s} \cdot \mu \cdot\left(1+\frac{\mu}{m}\right)+D_{n_{\mathrm{ss}}}+D_{f_{\mathrm{ss}}}\cdot \mu \cdot\left(1+\frac{\mu}{m}\right)\end{equation}$$ Under the assumption that degradation is null, this allows 8 to be simplified to the following: $$\begin{equation}P=f_{s s} \cdot \mu \cdot\left(1+\frac{\mu}{m}\right)\end{equation}$$

Calculating the Specific Growth Rate \(\mu\)

In order to calculate the specific growth rate, the following equation was applied to describe the binary fission of bacterial cells:$$\begin{equation}N=2^{\frac{t}{d_t}}\cdot N_0\end{equation}$$ Where \(N\) is population after time, \(N_0\) is the initial population size and t is time, and \(d_t\) is doubling time.

This can be rearranged to: $$\begin{equation}\log_{2}N=\frac{t}{t_d}+log{2}\cdot N_0\end{equation}$$ Given that: $$\begin{equation}\mu = \frac{(\ln{X} - \ln{X_0})}{t}\end{equation}$$ Where \(X_0\) is the initial cell concentration, \(X\) is the final cell concentration, while \(t\) is the time interval \(h\).

Then after one \(d_t\) , \(X\) = 2\(X_0\).

This allows rearrangement of the equation to:$$\begin{equation}\mu t_d = ln\:2\end{equation}$$Which was then applied to calculate the growth rates from the plot of \(log_{2} OD\) versus time to the experimental results.

To Improve

Within system modelling, assumptions had to be made about the models. These create limitations that shall be explored below.

Firstly, RNA polymerase and ribosomes are assumed to always be present in sufficient concentrations for the reaction to proceed, meaning there is no rate-limiting step to the reaction. If there were insufficient amounts of these elements, then this would slow down rates of transcription and translation meaning an extra step would need to be considered in future Simbiology models.

The process of protein binding and unbinding is also assumed to not be a time-limiting step in transcription and translation.

The rate-determining steps in the production of vsfGFP and the GS proteins are the transcription and translation of the genes being expressed. Other factors have been considered, such as the maturation rate of the GFP. The maturation rate of the GFP was important to consider as the folding of the GFP dictates the activation of its fluorophore.

That vsfGFP is so stable that during the 24 hour period of the experiment, no degradation occurs. If this assumption was to be changed, the Michaelis-Menten equation would be applied to calculate degradation rate.

The maturation constant was assumed from sfGFP, and this may differ due to the change in structure of the vsfGFP-0 used, particularly since vsfGFP is dimeric in contrast to the sfGFP monomer.

SimBiology Simulation and Compartment Analysis

In reference to Figure 1 of the structured modelling of vsfGFP expression, the simulation of the process was performed using Simbiology to estimate the amount of vsfGFP expressed from the induction of IPTG. The simulation was used as a visualisation technique to predict the overall performance within estimated kinetic parameters to understand how they affect the reaction process. Processes analysed include the uptake of IPTG across the cell membrane via active transport, to the translation of the mRNA that leads to the transcription of non-fluorescent vsfGFP \(GFP_{n} \). The \(GFP_{n} \) then matures to fluorescent GFP \(GFP_{f}\) that later makes the bacteria fluoresce. Due to the increase in bacterial growth, the concentration of vsfGFP also varies with time. Therefore using the Gompertz modelling equation as an estimate of bacterial growth, the amount of GFP is estimated.

In reference to Figure 3 below, the compartment with different species applied to the simulation of bacterial growth.

Figure 3: The growth compartment to simulate the growth of microbial cells that consists equationof the number of dormant bacteria which later divides at a specific activation rate. The model is used to estimate the cell concentration and at a certain activation rate, the number of dormant bacteria will either take a longer or a shorter time to undergo exponential growth.

The general equation used in computing the bacterial growth is given as; $$\begin{equation}\frac{dn_D}{dt} = -\alpha n_D\end{equation}$$ $$\begin{equation}\frac{dn_A}{dt} = r_a + r_g = \alpha n_D + \mu_{max}(1 - \frac{n_A}{N})n_A\end{equation}$$ Where \(\mu_{max}\) is the maximum growth rate (\(h^{-1}\)), \(\text{N}\) as the maximum concentration of the biomass at the point at which they are fed with limited nutrients, \(n_A\) as the concentration of active bacteria, \(n_D\) as the concentration of dormant bacteria, \(\alpha\) as the activation rate of the dormant cell (\(h^{-1}\)), and \(t\) is the time in hours \(\text{hr}\). The sum of populations of both the dormant bacteria and active bacteria \(n\), is then expressed as; $$\begin{equation}n = n_A + n_D\end{equation}$$ To finalise the simulation for the production of GFP, another compartment is added as illustrated in Figure 4 below.

Figure 4: The compartment of IPTG intake within a cell system. The reaction processes occurring within the compartment consist of the mass transfer of IPTG across the cell membrane (influenced by the co-transpport of lactose using the energy associated with the transmembrane proton motive force), the overall promoter activity, and the dilution and maturation rates of GFP proteins.

Due to the presence of the cell membranes and the lactose permease within of an E.coli bacteria, the mass transfer of IPTG is considered to occur due to the difference of the intercellular and extracellular difference of IPTG inducer concentration [4].Therefore, with the initial condition that the intercellular IPTG concentration is 0, the governing equation that can be used for IPTG uptake is described in Equation 17 below [4]: $$\begin{equation}\frac{\mathrm{d}\left [ IPTG \right ] }{\mathrm{d} t} = k_{c}\left (IPTG_{e} \right ) - \left (IPTG_{i} \right ) + K' \left [ IPTG_{e} \right ]\end{equation}$$ Where \([IPTG]_e\) is the extracellular IPTG concentration\(\mu\)M, \([IPTG]_f\) as the intracellular IPTG concentration \(\mu\)M, \(k_c\) as the mass transfer coefficient \(h^{-1}\), while \(K^{\prime} = \frac{k}{K_i}\) assuming that the IPTG initial concentration is less than \(K_i\).

Due to the increase in the number of active bacteria, it is predicted that the amount of fluorescent GFP increases with time as expressed in Equations below [7]; $$\begin{equation}\frac{d[GFP_n]}{dt} = gX_A(t) - k_m [GFP_n] -- \mu \cdot [GFP_n]- \frac{\gamma [GFP_n]}{[GFP_n] +[GFP_f] + M}\end{equation}$$ $$\begin{equation}\frac{d[GFP_f]}{dt} = gX_A(t) -k_m[GFP_f] -\mu \cdot [GFP_f]- \frac{\gamma [GFP_f]}{[GFP_n] +[GFP_f] + M}\end{equation}$$ Where \([GFP]_n\) is the concentration of non-fluorescent GFP, \([GFP]_f\) is the concentration of fluorescent GFP protein, \(M\) is the degradation capacity during the degradation process, \(g\) as the generation rate and \(m\) as the maturation constant and \(\gamma\) as the degradation rate [2].

The simulation was then performed with the assumptions that the degradation of GFP is negligible, the maturation constant to be the ln(2) divided by the maturation time of non-fluorescent GFP, which is assumed to be 13.6 min, the maximum concentration of biomass \(N\) as \(2\times 10^{9}\), the number of dormant bacteria to be 0.5 \times N, the mass transfer coefficient to be 0.213 \(h^{-1}\), \(K^{\prime}\) to be 0.0893 \(h^{-1}\)., \mu_max to be 1 while the generation rate becomes 0.76 \(h^{-1}\).

Results and Analysis

Following from Figure 5 below, the amount of \(GFP_f\) and \(GFP_n\) together with the biomass concentration is plotted with time.

Figure 5: The concentration profiles of the active and dormant bacteria with the amount of \(GFP_n\) and \(GFP_f\) where \(k_0\) is the maximum specific growth rate. \(\mu_{max}\).

It is observed that the number of the active bacteria first increases with time which then reduces. The amount of \(GFP_f\) and \(GFP_n\) increases over time until the maximum value. This is because of the influx of IPTG into the interstitial space of the bacteria at which the concentration becomes saturated.

Figure 6: The concentration profiles of the active and dormant bacteria with the amount of \(GFP_n\) and \(GFP_f\) where \(k_0\) is the maximum specific growth rate. \(\mu_{max}\). The maximum growth rate and the activation rate are manipulated to observe significant changes to the overall product formation.

When decreasing the activation rate, the time taken to produce GFP becomes greater compared to when the activation rate is 1.2. Shown in the Figure 6 above, the first run is when the activation rate constant becomes 0.98 while the 2nd run is when the activation rate becomes 1.2. When reducing the specific growth rate \(k_{0}\), it is expected that the slope of the bacterial culture growth becomes less steep. However, Iit is seen from the Figure above that the amount of active bacteria is more accumulated with lower growth rate. This could be that the value of \(\frac{\mu_{max}}{N}\) is < \(\alpha\) as shown in Figure 6 above.

This shows that the activation rate influences the time taken for the production of \(GFP_f\) where with an increasing activation rate, the less time \(GFP_f\) is producted. Also, the maximum growth rate has less effect to the growth of the microbial cells.
For the IPTG uptake, given that the initial concentration is 1 \(\mu\)M, it is assumed that the mass transfer coefficient is 0.213 and the overall profile for the IPTG intake becomes

Figure 7: The concentration profile of IPTG intake to the cell over time with change of the initial IPTG concentration.

Supporting Document of the Simbiology Simulation is available on iGEM Sheffield's GitLab

Biological Modelling References

[1] J. R. Kelly et al., “Measuring the activity of BioBrick promoters using an in vivo reference standard,” Journal of Biological Engineering, vol. 3, no. 1, p. 4, 2009, doi: https://doi.org/10.1186/1754-1611-3-4.
[2] J. H. J. Leveau and S. E. Lindow, “Predictive and Interpretive Simulation of Green Fluorescent Protein Expression in Reporter Bacteria,” Journal of Bacteriology, vol. 183, no. 23, pp. 6752–6762, Dec. 2001, doi: https://doi.org/10.1128/jb.183.23.6752-6762.2001.
[3] E. Frenzel, J. Legebeke, A. van Stralen, R. van Kranenburg, and O. P. Kuipers, “In vivo selection of sfGFP variants with improved and reliable functionality in industrially important thermophilic bacteria,” Biotechnology for Biofuels, vol. 11, no. 1, Jan. 2018, doi: https://doi.org/10.1186/s13068-017-1008-5.
[4] D. Calleja, A. Fernández-Castañé, M. Pasini, Carles de Mas, and Josep López‐Santín, “Quantitative modeling of inducer transport in fed-batch cultures of Escherichia coli,” Biochemical Engineering Journal, vol. 91, pp. 210–219, Oct. 2014, doi: https://doi.org/10.1016/j.bej.2014.08.017.
[5] [1]H. Alper, C. Fischer, E. Nevoigt, and G. Stephanopoulos, “Tuning genetic control through promoter engineering,” Proceedings of the National Academy of Sciences, vol. 102, no. 36, pp. 12678–12683, Aug. 2005, doi: https://doi.org/10.1073/pnas.0504604102.