Overview

---

In order to better understand the circuitry constructed for wet experiments and to some extent predict the performance of our system, providing feedback and guidance for the wet lab work, we have established a series of models, including a toxic gene expression model and a Cell Death Simulation Model.

The toxic gene expression model elucidates the principles of our system at the microscale, including the Epigenetic modification processes and the accumulation of toxic proteins over time. The Cell Death Simulation Model, on the macroscopic level, describes the process by which cells on the test paper come into contact with external solutions and alleviate inducer inhibition to initiate cell death. Both models are closely related to wet experiments, as illustrated in the diagram below:

Figure 1. The relationship between wet experiments and modeling

Toxic Gene Expression Model

---

In the toxic gene expression model , we employed a series of ordinary differential equations and Michaelis-menten equation to describe the wet experimental process. Subsequently, we utilized a solver to obtain numerical solutions for these systems of ordinary differential equations, deriving crucial numerical data from the wet experiments. Given that the wet experiments involved the study of yeast and Escherichia coli as host cells, we developed two similar yet distinct models based on the specific procedures and outcomes of the wet experiments.

Escherichia coli as the host

We introduced methyltransferase into Escherichia coli and induced its expression using IPTG. This enzyme performs methylation modifications on gene promoters, suppressing the expression of toxic genes and regulating cell death. To study the influence of IPTG concentration on methyltransferase levels, we developed a model for methyltransferase-induced expression.

Here, R represents the transcription product of the methyltransferase gene, while Met signifies its translation product. We employed a Hill function, denoted as f(S), to describe the influence of IPTG concentration, where S represents IPTG concentration, n is the Hill coefficient, and K is the half-maximal activation concentration. Additionally, we considered the production rate a and degradation rate b associated with the transcription product of the methyltransferase gene, as well as the production rate c and degradation rate d linked to its translation product.

When selecting the Hill function as the model, it's crucial to pay close attention and continually adjust the Hill coefficient n and half-maximal activation concentration K to ensure alignment with experimental data. This is because, when describing IPTG-induced gene expression, the choice of the Hill coefficient n influences the steepness of the gene expression response. A larger Hill coefficient indicates a smaller variation in gene expression response at lower IPTG concentrations, but a rapid increase in response at higher IPTG concentrations. Conversely, a smaller Hill coefficient suggests a relatively uniform gene expression response across the entire concentration range. Therefore, the choice of the Hill coefficient directly impacts the fit between the model and experimental data.

At the microscopic level, the process of methyltransferase-modified promoter binding to enzymes and substrates exhibits a remarkably high similarity. Inspired by this observation, we opted to employ the Michaelis-Menten equation to establish a model, specifically delineating the process of methyltransferase-modified promoter interaction.

Figure 2. The methyltransferase binds with the promoter.

In our research, Pro represents the unmodified promoter, while Met signifies the methyltransferase. Pro(Met) denotes the promoter that has undergone methylation by the methyltransferase. We utilize k_on to represent the rate of promoter methylation by the methyltransferase, and k_off represents the rate at which the modified promoter undergoes demethylation. [Pro]_T represents the total promoter concentration, which encompasses both unmodified and modified promoters.

For the sake of enhanced data processing and model analysis, we set the value of [Pro]_T to 1, while constraining the values of Pro within the range of 0 to 1. This restriction illustrates the proportion of unmethylated promoters to total promoters, represented as a percentage.

In order to delve deeper into the correlation between cell survival and death, we have taken into account not only the methylation levels but also essential factors such as toxic gene expression and the cell cycle. Our objective is to elucidate the nuanced balance inherent in these processes. Presented below are detailed descriptions of the toxic gene expression model and the cell growth model.

In this context, k2 represents the generation rate associated with transcription and translation processes, while the degradation rate of the product is denoted as k3. G represents the toxic gene expression product.

Here, e signifies the rate of cell self-proliferation, f represents the rate associated with cell apoptosis, g represents the impact factor of toxic gene G on cell survival, and C represents the cell count.

Results

Ultimately, we established interaction models depicting the variations over time in methyltransferase concentration (Met) and unmethylated modified promoter concentration (Pro). Additionally, we developed interaction models illustrating the protein levels of toxic gene expression and changes in cell count over time. These models were implemented using computer modeling tools. For specific details, please refer to the program descriptions in the "Software" section.

Figure 3. The trend of Met and Pro concentration over time. E.coli

Figure 4. The trend of toxic gene expression levels and cell numbers over time. E.coli

Through the analysis of the aforementioned two charts, it is evident that with the increase in methyltransferase expression levels, the proportion of unmodified promoters gradually decreases. The corresponding curve exhibits a rapid initial decline, eventually stabilizing. Simultaneously, at the initial stages, the expression of toxic genes synchronously increases with cell growth. However, as the expression of toxic genes reaches higher levels, the inhibitory effect on cell growth significantly intensifies. This is manifested in the graph by the occurrence of an inflection point in the cell growth curve. Subsequently, the number of cells undergoing growth gradually decreases, and the expression of toxic genes also stabilizes.

Yeast cells as the host

In our study focusing on yeast as the host cell, we employed a strategy involving targeted chromatin control through dCas9+ deacetylase to regulate the expression level of toxic genes. It is important to note that deacetylation acts on the three-dimensional chromatin structure, unlike methylation that targets the promoter region. Although deacetylation affects the overall gene structure, not limited to the promoter region, its mechanism of influencing gene expression can still be conceptually likened to a similar pathway. Therefore, we can approximate this situation by utilizing the previously described model used for Escherichia coli.

However, what sets this apart from the Escherichia coli model is that in our yeast experiments, the expression of deacetylase is constitutively active, and toxic genes are also in a constitutively active state. This is distinct from the situation in the Escherichia coli experiments where promoter methylation was considered.

Guidance on the results of wet experiments

When constructing the yeast cell growth model, Considering that the expression of deacetylase is constitutive, there was no need for an inducible expression model.Although gene expression in eukaryotic cells exhibits a more intricate complexity, when viewed from a simplified perspective, we have employed the model of methyltransferase-modified promoters in Escherichia coli and applied it to the context of deacetylase-modified chromatin.

In the actual experiments, we observed almost no colony growth on the culture medium initially. It was challenging to determine whether the insufficient effect was due to weak epigenetic modifications, the overpowering nature of toxic genes, or other factors. However, later, by employing the HML sequence recruitment method for epigenetic modifications, we observed colony growth on the culture medium. Hence, we inferred that the insufficient strength of epigenetic modifications was the cause.

To simulate this situation, we controlled the differential strength of epigenetic modification methods by adjusting the value of k1. Subsequently, we compared the simulation results with the wet lab experimental data and conducted further iterations. This process guided decision-making in the laboratory concerning the selection of epigenetic modification methods and the adjustment of parameters in dry experiments. Ultimately, we established an interaction model depicting the variation over time in protein content of toxic gene expression and cell count. The graphical representation of this model accurately reflected the real-world scenario.

Figure 5. The trend of toxic gene expression levels and cell numbers over time. Yeast (The solid line graph represents the results obtained when we applied dCas9+ deacetylase for epigenetic modification, whereas the dashed line graph corresponds to the data obtained using the HML sequence recruitment method for epigenetic modification.)

Further development from wet lab experiment

In our wet lab experiments, we observed a pronounced effect on the epigenetic modifications of weak promoters, while the modifications on strong promoters were relatively limited. During the initial modeling phase of the dry experiments, we initially assumed no significant correlation between k1 and k2. However, based on the results from our wet lab experiments, we identified a proportional relationship between k1 and k2, represented as k1=α*k2, where α signifies the degree of correlation between k1 and k2. Subsequently, we adjusted the model of the dry experiments by introducing the parameter α, replacing the previous control for k1, to precisely reflect the relationship between promoter strength and epigenetic modifications. We iteratively varied the value of α and compared the resulting graphs with the wet lab experimental outcomes, thereby determining the range of α values. Through this method, we delved deeply into the correlation between promoter activity and epigenetic modifications, providing detailed and accurate feedback for the wet lab experiments. The specific schematic diagram is presented below for reference.

Figure 6. The trend of Met and Pro concentration over time. Yeast (The solid line represents the graph generated when k2=4.236, while the dashed line corresponds to the graph obtained with k2=3.571.)

Figure 7. The trend of toxic gene expression levels and cell numbers over time. Yeast (The solid line represents the graph generated when k2=4.236, while the dashed line corresponds to the graph obtained with k2=3.571.)

Upon comparing the two charts, it is evident that as the promoter activity decreases (from solid line to dashed line), the epigenetic modification effect significantly intensifies. The percentage content of unmodified promoters decreases, leading to a substantial inhibition of toxic gene expression and a corresponding increase in cell survival rate. When k=3.571 (dashed line), toxic genes are no longer effectively suppressing cell growth; the entire cell growth curve shows a continuous upward trend without any inflection point (a different trend compared to the graph corresponding to k=4.236). With a mere 15% reduction in promoter activity, the proportion of unmodified promoters decreases by approximately 12%. Despite a mere 10% reduction in toxic gene expression, cell survival rate increases by almost 50%. This demonstrates the significant impact of toxic gene promoter activity on epigenetic modification effects, thereby eliciting pronounced effects on toxic gene expression and cell growth.

Cell Death Simulation Model

---

In order to simulate the macroscopic process of our test paper entering the environment to be tested, which involves the death of cells expressing toxic proteins, we have developed a Cell Death Simulation Model to describe this process. We simplified and improved upon Gang Cheng et al.'s tissue culture model, constructing a hybrid discrete-continuous (HDC) model suitable for our system. The HDC model consists of three components: cell distribution cellular automaton model, inducer diffusion model, and association model.

Cell Distribution Cellular Automaton Model:

We employ a cellular automaton model to simulate the process of cell proliferation and death on the test paper after the initiation of the suicide pathway in response to the inducer signal. Since different cells occupy different spatial positions within the hydrogel, and their exposure to the inducer in the external environment varies over time, it is necessary to model the distribution of cells within the hydrogel.

Here, we first need to address the initial distribution of cells in the hydrogel at t=0. Depending on how cells are introduced into the hydrogel, the initial distribution of cells can be categorized into two types: uniform distribution and surface distribution. Based on the steps involved in our hydrogel preparation, we have chosen to describe the initial distribution of cells using a uniform distribution.

The initial value of cell density is determined by the OD600 value of the cells and the steps involved in the preparation of the test paper, resulting in an initial cell density value of ρ₀ = 2.078x10⁸ cfu/μm3. Subsequently, the cell density ρ₀ is converted into the probability of cell appearance, p₀, at a certain location on the test paper using the formula p₀=(ρ₀ Vc)/V, where Vc is the average volume of cells, and V is the volume of the total space available for cell growth on the test paper. We use p₀ to randomly generate cell colonies in space and create an initial distribution map of cells in three-dimensional space:

Figure 8. Simulation of initial cell distribution on a hydrogel. Volume of (100x100x100) μm³

Inducer diffusion model:

During paper-based testing, the external solution only contacts the upper surface of the test paper. We take into account the diffusion of the inducer into the external solution on the test paper, which leads to a decrease in the inducer's concentration near the cells, initiating the expression of the suicide pathway. This process results in variations in the time at which cells at different locations on the test paper deactivate the inducer signal, causing different cells to activate the suicide pathway at different times.

To simplify calculations, we make the following assumptions in our model:

1. The inducer content in the external solution can be neglected. Since we are using IPTG as the inducer, which is predominantly present in the laboratory environment, we assume that the external environment being tested does not contain the inducer.

2. There is an ample amount of external solution, and during the diffusion process, the inducer concentration in the external solution can be considered as consistently 0.

3. The influence of convection on inducer diffusion is ignored. We assume that the temperature inside the hydrogel is uniform, and the inducer diffuses solely in the form of molecular diffusion.

4. We assume a constant temperature of T=298.15K (room temperature).

Without considering the effect of cells on the inducer on the test paper, the diffusion model follows the partial differential equation as follows:

Among these, C represents the concentration of the inducer at a specific location on the test paper, t is the time during which diffusion occurs, x, y, and z are three-dimensional spatial vectors, and De is the diffusion coefficient of the inducer. The value of De can be determined using an empirical formula:

Among these, Ф is the association factor of the solvent, M is the relative molecular mass of the solvent, T is the absolute temperature, μ is the viscosity of the solvent, and V_A is the molecular volume of the diffusing species. In our system, the solvent is hydrogel, so Ф is assigned a value of 2.6. T is set to room temperature, T=298.15K.

By substituting the calculated diffusion coefficient De into the partial differential equation, we obtain the diffusion results of the inducer over time. To facilitate visualization, we use a two-dimensional image to illustrate the diffusion process of the inducer on the vertical interface:

Figure 9. Inducer diffusion simulation. Where the concentration of the inducer is expressed by color.The closer the color is to red, the higher the inducer concentration.

Association model:

Next, we need to model the macroscopic process of the suicide pathway expression in different cells after deactivating the inducer signal. Whether the cell's suicide pathway can be successfully expressed depends on the concentration of inducer molecules around that cell. The inducer concentration, in turn, depends on spatial location and diffusion time. Therefore, it is necessary to integrate the first two models to construct an association model.

In this association model, we will establish a relationship between the inducer concentration, cell locations, and diffusion times. This will allow us to predict when and where the suicide pathway will be activated in different cells based on the inducer concentration profiles and spatial information.

The cell division rate r_s is calculated as follows:

Where r_s,max represents the cell's intrinsic maximum growth rate, Cn is the nutrient concentration at the cell's location, and K is the environmental carrying capacity for that cell.

The cell death rate r_d is determined by the results dC/dt of the toxic gene expression model.

Considering the limited computational power of our computer, we have made some simplifications to the model:

1. We set a threshold value C_i0, for when the suicide pathway starts to express and cells gradually die. In the calculation of the death rate r_d, C always takes the value C_i0.

2. We use the probability of new cells being generated at vacant positions ps instead of explicitly modeling the cell division rate r_s. We also use the probability of existing cells dying p_d instead of the cell death rate r_d. This significantly reduces our computational costs while still reflecting differences between cells to some extent.

Below is the complete iteration process diagram:

Figure 10. Flowchar of cellular automaton model iteration.

Finally, we conducted 200 iterations of the entire process, resulting in a simulation of the physiological changes of cells on the test paper after detecting the inducer signal.

Figure 11. The process of changes in the cells on the test paper after exposure to an external solution (Thickness = 30, C₀ = 1.0). Inducible substances enter the external solution from the upper surface of the test paper, causing the cells to express the suicide pathway sequentially from top to bottom, ultimately leading to cell death.

Feedback on wet experiments

Based on the simulation results from the cell death simulation model , we have observed that our whole-cell test paper, prepared according to our current experimental procedure, exhibits significant differences in the timing of cell death at different locations when used. Cells closer to the surface of the test paper quickly activate the suicide pathway, while cells near the bottom experience a slower decrease in the concentration of the inducer, resulting in a delayed expression of the suicide pathway. Consequently, these cells struggle to undergo cell death within the simulated timeframe. Based on the cell death simulation model , we speculate that it is possible to synchronize the activation of the suicide pathway in cells by reducing the thickness of the test paper or decreasing the initial inducer concentration, C₀, on the test paper.

Figure 12. The process of changes in the cells on the test paper after exposure to an external solution (Thickness =15, C₀=1.0).

Figure 13. The process of changes in the cells on the test paper after exposure to an external solution (Thickness =30, C₀=0.8).

The predictive results of the model indicate that by reducing the thickness of the test paper or decreasing the initial inducer content on the test paper, cells in the bottom layer of the paper can also undergo cell death within the simulated timeframe. This result provides an improvement strategy for the production steps of wet-lab whole-cell test paper.

How the model results affected the project design and development

---

On one hand, we have established a toxicity gene expression model that explains the microscopic processes of the induction and expression pathways of toxic genes. On the other hand, we have developed a cell death simulation model that describes the macroscopic process of cell deactivation of the inducer inhibition and initiation of cell death after the test paper comes into contact with an external solution. Our models not only provide a mechanistic understanding of how our suicide system operates but also offer valuable guidance for wet lab experiments.

The cell death simulation model suggests that for wet experiments, making the activation of the suicide pathway on the test paper more synchronized can be achieved by reducing the thickness of the test paper. In the toxicity gene expression model, we have provided control knobs for adjusting the values of each parameter (please refer to Software for more details). When switching between different strengths of promoters or epigenetic modification tools in wet experiments, you can conveniently adjust the corresponding values for predictions. Other teams using epigenetic modifications for gene expression regulation can also easily use our model for predictions.