Skip to content

Model

Section 1: Modelling N2O concentration

I. N2O dissolution

In order to predict N2O diffusion rates into the bacteria, first its expiration and dissolution rates must be modelled to establish the concentration gradient.

Expiration rate

Based on a paper by Dr. Hendrickx et al. at the Onze Lieve Vrouw Hospital, University of Stanford and Northwestern University (1), the composition of air expired by patients actively being administered Entonox as anaesthetic (~50% N2O) (0) was measured and plotted over time, showing this trend:

A graph showing the composition of air expired by patients actively being administered Entonox as anaesthetic over time. Graph taken from paper (1).

For the respective group, the graph shows that after the respective time of administering 4L/min of N2O, this was decreased to 0.4 L/min, and yet the expired N2O concentration decreased very slightly initially, then increased up to a potential maximum of ~60%. This suggests that the percentage concentration of expired N2O, stays roughly constant, and so can be modelled as such.

Dissolution rate

The dissolved concentration of N2O can be found with a given constant H_s^{cp} and partial pressure P by Henry’s Law:

C=H_s^{cp}P\\ P=c_{\%\N_2O}\times\ P_T

Here P_T is the total pressure, which for the canisters we used was 101.25kPa, and c_{%\N_2O} is the percentage concentration of N2O - 60%, as discussed above.

Temperature dependence of N2O solubility

H_s^{cp} was found to be 2.4\times{10}^{-4}\frac{mol}{m^3Pa} (2) at 298.15K, but to convert this to the hall constant at target temperature, H(T), the Van ‘t Hoff Equation is used:

\frac{d\left(\ln{H}\right)}{d\left(\frac{1}{T}\right)}=-\frac{\Delta_{sol}H}{R}

where T is the absolute temperature (K), and the right-hand side -\frac{\Delta_{sol}H}{R} is a constant that has been measured to be 2600K (2), so rearranging Eq. (1.3 ), letting -\frac{\Delta_{sol}H}{R}=k, where T_R and H_R are the respective temperature and hall constant at room temperature (298.15K):

\therefore\int_{H_R}^{H}{d(\ln{H})}=k\int_{T_R}^{T}d\left(\frac{1}{T}\right)\\\ \Longrightarrow\ln{\left(H\right)}-\ln{\left(H_R\right)}=\ln{\left(\frac{H}{H_R}\right)}=k\left[\frac{1}{T}-\frac{1}{T_R}\right]\\ \Rightarrow H(T)=H_R\exp{\left[k(\frac{1}{T}-\frac{1}{T_R})\right]}

Plugging this back into Henry’s Law:

C=H_R\exp{\left[k(\frac{1}{T}-\frac{1}{T_R})\right]}\times c_{\%\N2O}\times P_T\\ \therefore C=2.4\times{10}^{-4}\times\exp{\left[2600\left(\frac{1}{310}-\frac{1}{298.15}\right)\right]}\times60\%\times101250\\ \Longrightarrow C\approx10.45\frac{mol}{m^3}

II. N2O diffusion

When modelling our system for N2O breakdown, considering the diffusion of aqueous N2O into the E. coli cells is crucial - in fact, diffusion rates were a glaring issue with the first cycle of design.

Firstly, N2O crosses the bacteria surface membranes via simple diffusion, with seemingly no reliance on channel proteins (3) – unlike aquaporins for CO2 transport (4). Thus, this model will assume roughly constant diffusion across the phospholipid membrane, and only the outer one will be considered in this simplified model. So, we start with Fick’s Model of diffusion:

J=-D\frac{d\varphi}{dx}

Where J is the flux of the diffusing molecule (m^{-2}), D is the diffusion coefficient of the molecule (m^2s^{-1}) and \frac{d\varphi}{dt} is the difference in concentrations between ‘phases’.

This can be rewritten as:

J=-P_m\Delta c_s

Where P_m=\frac{K_pD}{d\ }, and thus:

J=\ -\frac{K_pD}{d}\Delta c_s (36)

Where K_p and d are the partition coefficient and thickness of the membrane respectively.

  • The diffusion coefficient D for N2O in water at 310K (the optimum temperature for nitrous oxide reductase (N2OR)) \approx2.5\times{10}^{-9}m^2s^{-1}(5)

  • The partition coefficient K_p could be found for only CO2 between water and membranes. The values were 0.95 at 300K (6), which we then estimate to be ~0.98 at 310K for N2O (see point 1 below), given the following assumptions:

    1. The temperature-dependence of the partition coefficient seems to follow a linear relationship, as K_p\propto T, as suggested in certain papers (11), (note this relationship was observed within the context of aerosols in air, but we extrapolate this to liquids and aqueous gas due to lack of literature surrounding these fluids. Thus, we can estimate that the K_p for CO2 at 310K would be \frac{310}{300}\times K_p^0\approx0.98
    2. The formula for the partition coefficient of CO2, K_p=\frac{\left[{CO}_2\right]\ in\ membrane}{\left[{CO}_2\right]\ in\ water} , is a ratio between the dissolved concentration in water and in the membrane. Since the solubility of CO2 and N2O in water is relatively similar (1.45g/L and 1.50 g/L respectively (7)) this ratio should be similar. Though this assumption ultimately reduces the accuracy of the final data, it is sufficient for a strong lower estimate of bacterial population.
    3. The partition coefficient was measured in egg lecithin liposome membranes, whose phospholipid composition is ~73% phosphatidylcholine (PC) (8), while E. coli outer membranes are ~20% phosphatidylglycerol (PG) and ~75% phosphatidylethanolamine (PE) (9). This may cause some effect on the value of the partition coefficients, however any variation found in the K_p for oxygen across different phospholipid membranes (DPPC, DMPC and DLPC) was within \pm10\times{10}^{-9} (10) and so will likely have a negligible effect on the final numerical result.
  • The thickness of the membrane d \approx 7.0 \times 10^{-9} , \text{m} (12)

  • The difference in solute concentration \Delta c_s will initially be \left[aqueous\ N_2O\right]-0=C from Eq. (1.9), as cytoplasm N2O concentration is assumed to be 0.

Thus, we can find the flux of N2O as:

\left|J\right|=0.98\times\frac{2.5\times{10}^{-9}}{7.0\times{10}^{-9}}\times10.45\approx3.7\frac{mol}{m^2}

We calculate the cross-sectional area A (assuming N2O is coming from one pointed direction, which we approximate to one end of a cylinder) of the E. coli membrane by modelling it as a cylinder. We take radius value of an E. coli bacterium to be 45\times{10}^{-6}\ m^2 (13).

A=\pi r^2\approx\pi\times\left(45\times{10}^{-6}\right)^2\approx6.4\times{10}^{-9}m^2

We now calculate number of N2O molecules (Q_{N_2O}) entering a single E. coli bacterium per second. We assume the N2O has dissolved homogeneously throughout the solution. This assumption is justified under our Cycle of Design 2 section.

Q_{N_2O}=J\times A\times N_A=3.7\times1.6\times{10}^{-14}\times6.022\times{10}^{23}\approx1.4\times{10}^{16}

Our simple model for the flow rate of small gas molecules transported from a gaseous to aqueous state, diffusing into a bacterial cell, as a function of temperature, can be summarised as:

Q_{\text{gas}}\left(T\right) = -\frac{K_p D}{d} \left(H_R \exp\left[k\left(\frac{1}{T} - \frac{1}{T_R}\right)\right] \cdot c_{\text{gas}} \cdot P_T \right) \cdot A \cdot N_A

Target population:

This now gives us an estimation of the maximum rate of diffusion into the cell, and so, given the maximum rate of reaction of the N2OR enzyme (v_{max}), we can estimate the minimum population size such that N2O would be constantly processed, which ensures maximum efficiency:

N_b=\frac{Q_{N2O}\left(mol\right)}{v_{max}}\\ =\frac{\frac{1.4\times{10}^{16}}{N_A}}{6.5\times{10}^{-20}}\approx3.6\times{10}^{11}bacteria

There are a few underlying assumptions to this model:

  1. The diffusion coefficient remains constant across water and phospholipid membranes, meaning we model the membrane as a slab of water except for the partition coefficient, which was determined experimentally across membranes and extracellular water (5).
  2. The diffusion coefficient of N2O is negligibly, or not at all, affected by pressure, as demonstrated experimentally in a paper from Imperial College London (Wang et al.) (5), making this seems a sensible assumption.
  3. The final result is clearly unrealistically large, though this is only in the instant at which N2O first dissolves into the solution, and the flux will change as the concentration in the cytoplasm increases and so the gradient into the cell decreases. Thus, we must consider the time-variation of the diffusion rate.

Time-variation of diffusion model

The initial Q_{N_2O} calculated must decrease as N2O begins to diffuse out of the cell as well as into it, so a time-dependent formula for \Delta c_s must be found. This can be done with Fick’s Second Law:

\frac{\partial\varphi}{\partial t}=D\frac{\partial^2\varphi}{{\partial x}^2}

In two dimensions, \frac{\partial^2\varphi}{{\partial x}^2} must be generalised to the Laplacian of concentration D\Delta\varphi:

\frac{\partial\varphi}{\partial t}=D\Delta\varphi

This has the general solution (14):

\varphi\left(x,t\right)=\frac{\alpha}{\sqrt t}\exp{\left[-\frac{x^2}{4Dt}\right]}

Where \alpha is a constant dependent on the boundary conditions, which we derive below. For a roughly ‘infinite’ N2O source with a constant concentration, this can be simplified:

Take an infinite bar with a thin, finite source to its left. The total number of solutes in the system must stay constant. Thus:

\int_{0}^{\infty}\varphi\left(x,t\right)dx=N

Where N is a constant that represents the total number of solutes. Thus, from Eq. (2.12):

\int_{0}^{\infty}{\frac{\alpha}{\sqrt t}\exp{\left[-\frac{x^2}{4Dt}\right]}dx}=N\\ =2\alpha\sqrt D\int_{0}^{\infty}\frac{1}{2\sqrt{Dt}}\exp{\left[-\left(\frac{x}{2\sqrt{Dt}}\right)^2\right]}dx\\ =2\alpha\sqrt D\int_{0}^{\infty}\exp{\left[-\left(\frac{x}{2\sqrt{Dt}}\right)^2\right]d\left(\frac{x}{2\sqrt{Dt}}\right)}\\

Let y=\frac{x}{2\sqrt{Dt}}. Then: N=2\alpha\sqrt D\int_{0}^{\infty}{e^{-y^2}dy\ }

Evaluating the Gaussian integral \int_{0}^{\infty}{e^{-y^2}dy\ }=\frac{\sqrt\pi}{2}

Leads us to: N=2\alpha\sqrt D\bullet\frac{\sqrt\pi}{2}\ \therefore\alpha=\frac{N}{\sqrt{\pi D}}

And so, plugging this back into Eq. (2.12): \varphi\left(x,t\right)=\frac{N}{\sqrt{\pi Dt}}\exp{\left[-\frac{x^2}{4Dt}\right]}

Now for a roughly infinite source with a constant concentration, we treat it as a bar of infinitely many of these thin sources of concentration \varphi_s and width du for some coordinate axis where -\infty<u<x. So, N=\varphi_sdu. Thus:

d(\varphi\left(x,t\right))=\frac{\varphi_sdu}{\sqrt{\pi Dt}}\exp{\left[-\frac{u^2}{4Dt}\right]}

To find the total diffusion source, we must sum all these thin sources:

\varphi\left(x,t\right)=\int_{x}^{\infty}{d(\varphi\left(x,t\right))}=\int_{x}^{\infty}{\frac{\varphi_s}{\sqrt{\pi Dt}}\exp{\left[-\frac{u^2}{4Dt}\right]}du\ }

Let y=\frac{u}{2\sqrt{Dt}}. Then:

\varphi\left(x,t\right)=\frac{\varphi_s}{\sqrt\pi}\int_{\frac{x}{2\sqrt{Dt}}}^{\infty}{\frac{2\sqrt{Dt}}{\sqrt{Dt}}e^{-y^2}dy}=\frac{2\varphi_s}{\sqrt\pi}\int_{\frac{x}{2\sqrt{Dt}}}^{\infty}{e^{-y^2}dy}\\ =\frac{2\varphi_s}{\sqrt\pi}\left\{\int_{0}^{\infty}{e^{-y^2}dy}-\int_{0}^{\frac{x}{2\sqrt{Dt}}}{e^{-y^2}dy}\right\}\\ =\frac{2\varphi_s}{\sqrt\pi}\left\{\frac{\sqrt\pi}{2}-\int_{0}^{\frac{x}{2\sqrt{Dt}}}{e^{-y^2}dy}\right\}=\varphi_s\left\{1-\frac{2}{\sqrt\pi}\int_{0}^{\frac{x}{2\sqrt{Dt}}}{e^{-y^2}dy}\right\}\\ =\varphi_s\left\{1-\text{erf}{\left[\frac{x}{2\sqrt{Dt}}\right]}\right\}\\ \therefore\ \varphi\left(x,t\right)=\varphi_s\left\{1-\text{erf}\left[\frac{x}{2\sqrt{Dt}}\right]\right\}

Where erf is the error function, which was only calculated numerically using Desmos for the graphs used. Plugging in the values for aqueous N2O concentration and the diffusion coefficient of our system:

\varphi\left(x,t\right)=10.45\times\left\{1-\text{erf}{\left[{10}^4\times\frac{x}{\sqrt t}\right]}\right\}

III. N2OR synthesis rate

Now having a profile for the N2O diffusion into the bacteria, we create a model for the rate of the nitrous oxide reductase (N2OR) reaction to inform design decisions on how to optimise the system.

Transcription rate & choosing a promoter

The transcription rate of N2OR must be determined to find the enzyme concentration in the transformed E. coli. The Mass Action Kinetics model is used for this: We decided on a constitutive promoter was as constant N2OR production is required for our N2O breakdown design to be used continuously. Thus, there are no transcription factors that need to be modelled, leading us to adopt this model from this paper (15):

B=\frac{B_{max}\times\left[RNAP\right]}{K_d+\left[RNAP\right]}

Where B is the proportion of DNA of the target gene that is bound by RNAP, B_{max} is the maximum binding for a given promoter, K_d (M) is the dissociation constant and \left[RNAP\right] is the concentration of RNA Polymerase (RNAP) in the cell. We took several considerations into account when choosing the promoter:

  1. Clearly, the transcription rate B, represented by the proportion of DNA bound to RNAP, is proportional to B_{max}. Thus, B_{max} must be maximized when choosing the promoter. Several of the promoters in the J231xx family had a B_{max} of 0.99, namely: J23119, J23108, and J23115 (15).

  2. UP susceptibility: UP (upstream) elements can affect transcriptions by altering the binding of RNAP to the DNA. Since UP elements would not be considered for the genetic design of this project, a promoter with high performance assuming there are no UP elements was desired. This was judged using the following graph (15):

A graph of the mean relative unit of GFP luminescence using several different promoter in the J231XX family. Taken from (15).

Thus:

B=\frac{0.99\times2.5\times{10}^{-6}}{37.78\times{10}^{-9}+2.5\times{10}^{-6}}\approx0.98

J23199 is the clear candidate for this criterion, as it is the most active promoter with and without upstream elements.

  1. K_d: the dissociation constant is also a factor in the transcription rate, but its effect is much less significant, as B\propto\frac{1}{K_d+\left[RNAP\right]} and so it will have a limited effect since K_d\ll[RNAP] (~2 orders of magnitude). Since K_d varied by ~2.5\times{10}^{-8}M (15), this is an even more insignificant consideration and can be ignored.

Therefore, J23119 was chosen as the promoter for the genetic design, having the highest B_{max} of 0.99 and activity measured by relative fluorescence. Thus, with its K_d being ~37.78\times{10}^{-9} (15):

B=\frac{0.99\times\left[RNAP\right]}{37.78\times{10}^{-9}+\left[RNAP\right]}

The DH5α strain of Escherichia coli, growing in carbon sources, has a doubling time of 45 minutes at 310K (16) and average cell volume of \approx(2.0\pm1.3)\times{10}^{-18}m^3, or 2.0 fL (17). We assume to be roughly equivalent to that of the E. coli strain B/rA which was measured to be 2.5\times{10}^{-6}M𝑀 at 310K, with a doubling time of 24 minutes (18). The justifications for this assumption are:

  1. A paper by Shepherd et al. from the Journal of Bacteriology (18) found that \left[RNAP\right] seemed independent of doubling time, and so this disparity between DH5α chosen for our project and B/rA used in this same paper, can be ignored.
  2. This B/rA strain has an average cell volume of 1.25\times{10}^{-18}m^3 (18), which falls within the first standard deviation of the value for DH5α, thus suggesting that the results measured in the paper could reasonably be taken from cells of the same volume as DH5α, and so this is not sufficiently significant to discount the validity of the assumption.

Thus:

B=\frac{0.99\times2.5\times{10}^{-6}}{37.78\times{10}^{-9}+2.5\times{10}^{-6}}\approx0.98

We can then use the Hill Equation, that for a reaction A+n\times B\rlhar C:

\frac{d\left[C\right]}{dt}=k_a\left[A\right]\left[B\right]^n-k_d[C]

And at equilibrium:

\frac{d\left[C\right]}{dt}=0\\ \Rightarrow k_a\left[A\right]\left[B\right]^n=k_d[C]\\ \Rightarrow K_d=\frac{k_d}{k_a}=\frac{\left[A\right]\left[B\right]^n}{\left[C\right]}

Therefore, the proportion of molecules A that are bounded (in our case, DNA), B is such that:

B=\frac{\left[C\right]}{\left[A\right]+\left[C\right]}\\ \therefore\left[C\right]=[mRNA]=A+C∙B=RNAP∙B

So, plugging in the values established above:

\left[mRNA\right]=2.5\times{10}^{-6}M\times0.98\approx2.4\times{10}^{-6}M Before the mRNA is translated, its natural decay rate must be considered. The median half-life of *E. coli* mRNA molecules was found to be 5 minutes (300s) (19), and thus representing half-life as \tau: \left[mRNA\right]^\prime=\left[mRNA\right]\times2^{-\frac{t}{\tau}}\\ \Rightarrow\frac{d\left[mRNA\right]\prime}{dt}=-\frac{\ln{2}}{\tau}\left[mRNA\right]\times2^{-\frac{t}{\tau}}\\ =-\frac{\ln{2}\ }{\tau}\left[mRNA\right]\prime\approx-2.31\times{10}^{-3}[mRNA]

Having assumed equilibrium in the calculation for [mRNA], and thus no net increase in [mRNA], this subsequently seems to suggest total [mRNA] will decrease. However, factoring in the mRNA transcription rate, when equilibrium is disrupted, renders this a negligible decrease compared to total [mRNA], and thus it being ignored.

Translation rate We can represent the translation process as mRNA\rightarrow mRNA+Protein, because mRNA molecules can be reused after translation by a ribosome.

Then, during translation, the rate of protein synthesis can be written as:

\frac{d\left[Protein\right]}{dt}=k_{TL}[mRNA]

Where k_{TL} isthe translation rate constant. We could not find the value of this constant despite literature reviews, leading us to model its value. For this calculation, we used the equation (20):

Where:

TA=\frac{\left[mRNA\right]\bullet r i b d e n\bullet r i b o c c}{K_m+\left[mRNA\right]\bullet r i b d e n\bullet r i b o c c}
  • ribden is the ribosome density along an mRNA molecule and is usually found empirically. Unfortunately, we could not find data regarding NosZ but found the ribosome density for to be 0.64 (20). The accuracy of this value for NosZ cannot be evaluated, and must be treated at most as a ‘ballpark’ figure. Indeed, ribosome density cannot be greater than 1, and the median in yeast is 0.36 and can be considered as the lower bound (21), so a sensible uncertainty for this value is posited as 0.64 ± 0.28.
  • ribocc is the ribosome occupancy, the percentage of ribosomes actively translating the mRNA- which on average in E. coli in stationary states was found to be 91% ± 7% (21). A lower bound of 84% will be used to provide a minimum rate of reaction, and upper bound 98% for a maximum rate. This will hopefully tackle the potential inaccuracy of using the value averaged over the entire E. coli genome rather than for NosZ specifically.
  • The K_m has been estimated previously by maximising the Spearman’s Rank Correlation Coefficient with experimental data of N2OR concentrations with the predicted value using this calculation (20). The maximum correlation was found when K_m=0.06, for the E.coli genome.

k_p and the translation elongation rate (\eta), which is the rate at which a polypeptide is ‘elongated’ by th e addition of another amino acid to the chain during translation, are equivalent (21). As such, X was used in place of k_p.

k_{TL}=\eta\bullet\frac{\left[mRNA\right]\bullet r i b d e n\bullet r i b o c c}{K_m+\left[mRNA\right]\bullet r i b d e n\bullet r i b o c c}

For E. coli growing at a rate of <0.035\ h^{-1}, \eta\approx8aa\ s^{-1} (22). In an ideal, equilibrated bacterial population, the net growth rate is 0\ h^{-1}. Though this will not be true for each individual bacterium, we use this as the limit at very small growth rates to simplify our model. This comes from these graphs, taken from (23).

A graph of translation elongation rate (TEA) (aa/s) against growth rate (/h). The curve shows a projection of the translation elongation rate at a growth rate of exactly 0/h: 8aa/s. It can also be seen that at very low variation in growth, the effect on the TEA is somewhat negligible. Paper taken from (20).

The NosZ gene which codes for N2OR contains 151 amino acids, and thus the time taken to translate one by one ribosome is:

T=\frac{151}{9}\approx16.8s

We add to this the initiation time - the time taken for a ribosome to successfully attach to the mRNA molecule. It is quoted as ~10s for E. coli and is independent of growth rate and temperature, rendering it likely very accurate (22):

T=16.8+10=26.8s

Thus, per ribosome, the rate of translation is:

\frac{1}{T}\approx0.037s^{-1}

Plugging these values into Eq. (3.17), using [mRNA] calculated in Eq. (3.10) the lower and upper bounds are:

k_{TL\ LB}=0.037\times\frac{\left[2.4\times{10}^{-3}\right]\times0.36\times0.84}{0.06+\left[2.4\times{10}^{-3}\right]\times0.36\times0.84}\\ \approx4.4\times{10}^{-4}s^{-1}\\ k_{TL\ UB}=0.037\times\frac{\left[2.4\times{10}^{-3}\right]\times0.92\times0.98}{0.06+\left[2.4\times{10}^{-3}\right]\times0.92\times0.98}\\ \approx1.3\times{10}^{-3}s^{-1}

The upper bound for the rate is ≈290% greater than the lower bound, showing significant uncertainty and thus a need for a more accurate model, or empirical data from direct study of the NosZ ge ne. We use these bounds further down our model.

Proteolysis

Proteins eventually undergo proteolysis, where they are broken down and/or inactivated. This provides a dissociation reaction for equilibrium such that:

\frac{d\left[Protein\right]}{dt}=k_{TL}\left[mRNA\right]-k_X[Protein]

Where k_X is the rate of proteolysis. And thus, in equilibrium we theorised that:

\frac{d\left[Protein\right]}{dt}=0\\ \Longrightarrow k_{TL}\left[mRNA\right]=k_X[Protein]\\ \Longrightarrow\left[Protein\right]=\frac{k_{TL}}{k_X}\left[mRNA\right]

For N2OR, proteolysis (inactivation) rate was experimentally estimated as 0.4h^{-1}\approx1.1\times{10}^{-4}s^{-1} (25), in aerobic conditions and thus in the presence of potentially non-competitive inhibitors such as azides and oxygen. As such, roughly the same conditions for N2OR inactivation exist in our aerobic system.

Plugging this into Eq. (3.14), the lower bound is:

\left[N_2OR\right]=\frac{4.4\times{10}^{-4}}{1.1\times{10}^{-4}}\times2.4\times{10}^{-3}\approx9.6\times{10}^{-3}\frac{mol}{m^3}

And the upper bound:

\left[N_2OR\right]=\frac{1.3\times{10}^{-3}}{1.1\times{10}^{-4}}\times2.4\times{10}^{-3}\approx2.8\times{10}^{-2}\frac{mol}{m^3}

The model we developed to find the protein concentration in a stable bacterial cell can be thus summarised as:

\left[Protein\right]=\frac{1}{k_{X\ \ }}\left(\eta\bullet\frac{\left[mRNA\right]\bullet r i b d e n\bullet r i b o c c}{K_m+\left[mRNA\right]\bullet r i b d e n\bullet r i b o c c}\right)\left(\frac{B_{max}\times\left[RNAP\right]^2}{K_d+\left[RNAP\right]}\right)

IV. N2OR rate of reaction

Given the concentration of the N2OR enzyme, we model rate of reaction using the Michaelis-Menten equation:

For a reaction of the form:

E+S\rlhar ES\rightarrow E+P

Where E (enzyme) is N2OR, S (substrate) is N2O, ES is the enzyme-substrate complex, and P (product) is N2 and H2O, we use the Michaelis-Menten equation:

\frac{d\left[P\right]}{dt}=k_{cat}\left[E_0\right]\frac{\left[S\right]}{K_d+\left[S\right]}

k_{cat}is the catalysis constant, and \left[E_0\right] is the total concentration of enzymes, bound and unbound, and k_{cat}\left[E_0\right] is also called the maximum rate of reaction v_{max}. Thus, this can equivalently be written as:

v=\frac{v_{max}\left[S\right]}{K_m+\left[S\right]}

Note, however, this would be an extreme simplification of the N2O breakdown reaction, as we have ignored the necessity of 2 electrons and 2 protons needed for the full reaction:

N_2O+2e^-+2H^+\rightarrow H_2O+N_2

These can be treated as another 2 substrates in the reaction, e- and H^+, and so the Michaelis-Menten equation (Eq. (4.3)) can be extended to a 3-substrate form (26):

v=\frac{\left[E_0\right]}{\frac{1}{k_0}+\frac{1}{k_A\left[A\right]}+\frac{1}{k_B\left[B\right]}+\frac{1}{k_C\left[C\right]}+\frac{1}{k_{AB}\left[A\right]\left[B\right]}+\frac{1}{k_{AC}\left[A\right]\left[C\right]}+\frac{1}{k_{BC}\left[B\right]\left[C\right]}+\frac{1}{k_{ABC}\left[A\right]\left[B\right]\left[C\right]}}

Where \left[E_0\right] is the enzyme concentration and k_0 is a rate constant such that k_0\left[E_0\right]=v_{max}. Thus, we simplify the above to:

v=\frac{v_{max}\left[A\right]\left[B\right]\left[C\right]}{\left(K_m^a+\left[A\right]\right)\left(K_m^b+\left[B\right]\right)\left(K_m^c+\left[C\right]\right)}

However, since there are two e^- and H^+ we must treat these as 2 separate reactants:

v=\frac{v_{max}\left[A\right]\left[B\right]^2\left[C\right]^2}{\left(K_m^a+\left[A\right]\right)\left(K_m^b+\left[B\right]\right)^2\left(K_m^c+\left[C\right]\right)^2}

Written explicitly:

v=\frac{v_{max}\left[N_2O\right]\left[H^+\right]^2\left[e^-\right]^2}{\left(K_m^{N2O}+\left[N_2O\right]\right)\left(K_m^{H+}+\left[H^+\right]\right)^2\left(K_m^{e-}+\left[e^-\right]\right)^2}

H^+

For H^+, the pH of the bacterial cytoplasm was used to estimate the average concentration of H^+ ions. The normal pH of E. coli ranges from 7.2 to 7.8 (27), and the optimum pH for N2OR can be estimated experimentally, as visible in this graph (28):

A graph of the different stages of denitrificaiton and their maximum reduction rate at different pH. It can be seen that the maximum rate of nitrous oxide breakdown is somewhere between pH 7.75-8.25. Taken from (28)

N2O reduction is maximised between pH 7.75-8.25 (approximating this to a smooth curve). Combined with the pH range found for E. coli, a pH of 7.8 seems like a sensible estimate for the optimum pH for the bacterial system to model. To find [H+]:

\left[H^+\right]={10}^{-pH}\\ \approx{10}^{-7.8}\approx1.6\times{10}^{-8}mol\ L^{-1}=1.6\times{10}^{-5}mol\ m^{-3}

As for K_m^{H+}, though no direct data can be found in the case of N2OR activity, other proton donors in enzymatic reactions have been tested. For Uroporphyrinogen decarboxylation, the proton dissociation constant was measured to be ~1.7\ \times\ {10}^{-24}M=1.7\times{10}^{-21}mol\ m^{-3} (29). Though this value seems dubiously small, it provides evidence that K_m^{H+}\ll[H^+] and thus K_m^{H+}+\left[H^+\right]\approx[H^+].

e-

For \left[e^-\right], the concentration of the predominant electron donors can be used as an estimate. The main electron donors in denitrification are theorised to be Fe2+ and Mn2+, accounting for 97% of the electron contribution (30). Thus, assuming each ion donates 2 electrons to form the 2+ ions, \left[e^-\right] can be estimated as 2\left[{Fe}^{2+}\right] +2\left[{Mn}^{2+}\right].

  • \left[{Fe}^{2+}\right]\approx1\mu M=1\times{10}^{-3}\ mol\ m^{-3} (31)
  • \left[{Mn}^{2+}\right]\approx5\mu M=5\times{10}^{-3}\ mol\ m^{-3} (32)

Thus \left[e^-\right]\approx2\times\left(1+5\right)\times{10}^{-3}=1.2\times{10}^{-2}\ mol\ m^{-3}

As for K_m^{e-}, we can only estimate this using data from other reactions, due to the lack of experimental data modelling the electron donations in N2O reduction. However, values for other electron donors to other enzymes can be used to give a range: in CYP2C9, enzymes which require electron donors such as NADPH, the K_m ranges from 33nM to 141nM=0.33\times{10}^{-6}mol\ m^{-3} (33). We consider the lower potential bound _2O degradation, to constrain the worst-case-scenario, and find the greatest possible difference between K_m and \left[H^+\right]. As such, we use 1.4\times{10}^{-6}mol\ m^{-3}.

N2O

For \left[N_2O\right], we will use the function derived in Section II. to provide a time-dependent rate of reaction, reflecting how it would decrease as N2O concentration did too. The x-value would be taken from the centre of a cell and would represent the average concentration throughout the cytoplasm. E. coli is on average 2.47\times{10}^{-6}m long (13), and so for \left[N_2O\right], \varphi(1.2\times{10}^{-6},t) will be used.

V_{max} and K_m

v_{max} was measured to be 1.73\pm\ 0.05\ \mu mol per mg of protein per minute \approx2.8\times{10}^{-5}mol{\ s}^{-1} per g of protein (34). This is equivalent to the actual v_{max} divided by the mass of enzyme available: k=2.8\times{10}^{-5}mol{\ s}^{-1}\ g^{-1} such that:

v_{max}=k\bullet\left[N_2OR\right]\bullet{vol}_{cell}\bullet m_r

Where m_r is the relative molecular mass \approx120kDa=1.2\times{10}^5g\ mol^{-1} (35), and thus, using the average cell volume of E. coli as 2.0\times{10}^{-18} (17), the upper bound for v_{max} is, using \left[N_2O\right] from Eq. (3.27):

v_{\max{UB}}=2.8\times{10}^{-5}\times2.8\times{10}^{-2}\times2.0\times{10}^{-18}\times1.2\times{10}^5\\ \approx1.9\times{10}^{-19}mol{\ s}^{-1}

And the lower bound, using \left[N_2O\right] from Eq. (3.28) is:

v_{\max{LB}}=2.8\times{10}^{-5}\times9.6\times{10}^{-3}\times2.0\times{10}^{-18}\times1.2\times{10}^5\\ \approx6.5\times{10}^{-20}mol\ s^{-1}

K_m was measured to be 25.6\pm2.3\ µM\approx2.6×10^{-2} mol\ m^{-3} (34).

Calculating rates:

Starting with the simplified single-substrate model:

v_1=\frac{v_{max}\left[N_2O\right]}{K_m+\left[N_2O\right]}

We can write the expanded model as a multiple of the above due to the assumed constant nature of the other substrate concentrations:

v_2=\frac{v_{max}\left[N_2O\right]\left[H^+\right]^2\left[e^-\right]^2}{\left(K_m^{N2O}+\left[N_2O\right]\right)\left(K_m^{H+}+\left[H^+\right]\right)^2\left(K_m^{e-}+\left[e^-\right]\right)^2}\\ =\mathcal{K}v_1

Where:

\mathcal{K}=\frac{{(\left[H^+\right]}^2\left[e^-\right]^2)}{\left(K_m^{H+}+\left[H^+\right]\right)^2\left(K_m^{e-}+\left[e^-\right]\right)^2}\\ \mathcal{K}=\frac{\left(1.6\times{10}^{-5}\right)^2\times\left(1.2\times{10}^{-2}\right)^2}{\left(1.6\times{10}^{-5}\right)^2\times\left(1.2\times{10}^{-2}+1.4\times{10}^{-6}\right)^2}\\ \approx0.99977

The expanded model therefore does suggest that the simplified one is an overestimate of ~0.023%. However, this is negligible in comparison with the 290% variation between the upper and lower bound for v_{max} and thus the rate of reaction itself.

Evaluating v_1

First the lower bound:

v_{1\ LB}\left(t\right)=6.5\times{10}^{-20}\times\frac{\varphi(1.2\times{10}^{-6},t)}{2.6\times{10}^{-2}+\varphi(1.2\times{10}^{-6},t)}

And:

v_{1\ UB}\left(t\right)=1.9\times{10}^{-19}\times\frac{\varphi(1.2\times{10}^{-6},t)}{2.6\times{10}^{-2}+\varphi(1.2\times{10}^{-6},t)}

Clearly, for \varphi\gg2.6\times{10}^{-2},\ v=v_{max}, and so this will be the maximum rate in our system from which we can evaluate minimum bacterial population (Section II.)

Limitations of this model:

  • Modelling of e^- and H^+ as ‘substrates’ of the enzyme likely leads to inaccuracies in our final result, as subtler complexities such as the oxidation of the N2OR copper sites have been ignored for simplicity. This could be a point for a future team to improve on our model.
  • Our estimation of \left[e^-\right] is inevitably an underestimate, as other ions and compounds will contribute electrons to the reaction, but 97% of electron contribution in denitrification was deemed a sufficiently accurate representation for this model, which already incorporates much greater uncertainty in other data values.
  • The data taken for these values was taken from different environments for the E. coli colonies, and so the contributions of different ions to denitrification may differ based on availability. We had to ignore this for the sake of our model due to a lack of alternative data.
  • The need to estimate the values for the dissociation constants K_m ultimately negates the benefit of the increased accuracy of the expanded model to the actual reaction. We would welcome future teams to investigate and attempt to experimentally ascertain these constants for improved models of N2O breakdown.

V. Combined model for [N2O]

Combining the rate of reaction with the rate of diffusion, one model can be formed:

\frac{d\left[N_2O\right]}{dt}=-N_bv\left[N_2O\right]+JA\\ =-N_b\mathcal{K}\frac{v_{max}\left[N_2O\right]^2}{K_m+\left[N_2O\right]}-D\frac{\partial\varphi}{\partial x}A

Where is N_bthe number of bacteria in the population. Since we have previously defined as \left[N_2O\right] as \varphi\left(x,t\right):

\frac{\partial}{\partial t}\left(\varphi\left(x,t\right)\right)=-N_b\mathcal{K}v_{max}\frac{(\varphi{\left(x,t\right))}^2}{K_m+\varphi\left(x,t\right)}-DA\frac{\partial}{\partial x}\left(\varphi\left(x,t\right)\right)\\ \therefore\varphi\left(x,t\right)=-N_b\mathcal{K}v_{max}\int_{0}^{T}{\frac{(\varphi{\left(x,t\right))}^2}{K_m+\varphi\left(x,t\right)}dt}-DA\int_{0}^{T}{\frac{\partial}{\partial x}\left(\varphi\left(x,t\right)\right)dt\ }

Firstly, \frac{\partial}{\partial x}\left(\varphi\left(x,t\right)\right) can be found:

\frac{d\varphi}{dx}=\frac{d}{dx}\left(\varphi_s-\frac{2\varphi_s}{\sqrt\pi}\int_{0}^{\frac{x}{2\sqrt{Dt}}}{e^{-y^2}dy}\right)

Letting u=\frac{x}{2\sqrt{Dt}}, such that \frac{du}{dx}=\frac{1}{2\sqrt{Dt}}:

=0-\frac{2\varphi_s}{\sqrt\pi}\frac{d}{du}\left(\int_{0}^{u}{e^{-y^2}dy}\right)\bullet\frac{du}{dx}\\ =-\frac{2\varphi_s}{\sqrt\pi}\bullet\frac{1}{2\sqrt{Dt}}\bullet\frac{d}{du}\left(\int_{0}^{u}{e^{-y^2}dy}\right)

By the Fundamental theorem of calculus: \int_{a}^{b}{f^\prime\left(t\right)}dt=f\left(b\right)-f\left(a\right) If e^{-y^2}=f^\prime\left(y\right) a=0,\ b=u :

\int_{0}^{u}{f^\prime\left(y\right)}dy=f\left(u\right)-f\left(0\right)\\ \Rightarrow\frac{d}{du}\int_{0}^{u}{f^\prime\left(y\right)}dy=\frac{d}{du}\left(f\left(u\right)-f\left(0\right)\right)\\\ =f^\prime\left(u\right)=e^{-u^2}\\\ \therefore\frac{d}{du}\left(\int_{0}^{u}{e^{-y^2}dy}\right)=e^{-u^2}

Plugging this back into Eq. (5.6):

\frac{d\varphi}{dx}=-\frac{2\varphi_s}{\sqrt\pi}\bullet\frac{1}{2\sqrt{Dt}}\bullet\exp{\left[-\left(\frac{x}{2\sqrt{Dt}}\right)^2\right]}\\ =-\frac{\varphi_s}{\sqrt{\pi Dt}}\exp{\left[-\frac{x^2}{4Dt}\right]}

So,

\int_{0}^{T}{\frac{\partial}{\partial x}\left(\varphi\left(x,t\right)\right)dt\ }=-\frac{\varphi_s}{\sqrt{\pi D}}\int_{0}^{T}{\frac{1}{\sqrt t}\exp{\left[-\frac{x^2}{4Dt}\right]}}dt

Letting u=\exp{\left[-\frac{x^2}{4Dt}\right]} and \frac{dv}{dt}=\frac{1}{\sqrt t} such that \frac{du}{dt}=\frac{x^2}{4Dt^2}\exp{\left[-\frac{x^2}{4Dt}\right]} and v=2\sqrt t, by parts:

\int_{0}^{T}{\frac{1}{\sqrt t}\exp{\left[-\frac{x^2}{4Dt}\right]}}dt=\left.2\sqrt t\exp{\left[-\frac{x^2}{4Dt}\right]}\right|_0^T-\int_{0}^{T}{\frac{x^2\exp{\left[-\frac{x^2}{4Dt}\right]}}{2Dt^\frac{3}{2}}dt}

Let w=\frac{x}{2\sqrt{Dt}} such that \frac{dw}{dt}=-\frac{x}{4\sqrt D\ t^\frac{3}{2}\ }

\int_{0}^{T}{\frac{x^2\exp{\left[-\frac{x^2}{4Dt}\right]}}{2Dt^\frac{3}{2}}dt}=\frac{2x}{\sqrt D}\int_{0}^{\frac{x}{2\sqrt{DT}}}{e^{-w^2}dw}\\ =\frac{x\sqrt\pi}{\sqrt D}\bullet\frac{2}{\sqrt\pi}\int_{0}^{\frac{x}{2\sqrt{DT}}}{e^{-w^2}dw}=\frac{x\sqrt\pi}{\sqrt D}\text{erf}{\left(\frac{x}{2\sqrt{DT}}\right)}

Plugging this back into Eq. (5.16): \int_{0}^{T}{\frac{1}{\sqrt t}\exp{\left[-\frac{x^2}{4Dt}\right]}}dt=\left.2\sqrt t\exp{\left[-\frac{x^2}{4Dt}\right]}\right|0^T-\frac{x\sqrt\pi}{\sqrt D}\text{erf}{\left(\frac{x}{2\sqrt{DT}}\right)}\ =2\sqrt T\exp{\left[-\frac{x^2}{4DT}\right]}-\frac{x\sqrt\pi}{\sqrt D}\text{erf}{\left(\frac{x}{2\sqrt{DT}}\right)}\ \int{0}^{T}{\frac{\partial}{\partial x}\left(\varphi\left(x,t\right)\right)dt\ }=\frac{\varphi_s}{\sqrt{\pi D}}\left{\frac{x\sqrt\pi}{\sqrt D}\text{erf}{\left(\frac{x}{2\sqrt{DT}}\right)}-2\sqrt T\exp{\left[-\frac{x^2}{4DT}\right]}\right}

Secondly,

\int_{0}^{T}{\frac{(\varphi{\left(x,t\right))}^2}{K_m+\varphi\left(x,t\right)}dt} cannot be analytically solved, as far as we tried, but some estimations can be made: For \varphi\left(x,t\right)\gg\ K_m, which is certainly true initially; 10.45\gg0.02:

\int_{0}^{T}{\frac{(\varphi{\left(x,t\right))}^2}{K_m+\varphi\left(x,t\right)}dt}\approx\int_{0}^{T}\varphi\left(x,t\right)dt\\ =\int_{0}^{T}{\varphi_s\left\{1-\text{erf}{\left[\frac{x}{2\sqrt{Dt}}\right]}\right\}dt}=\varphi_sT-\varphi_s\int_{0}^{T}{\etext{erf}rf{\left[\frac{x}{2\sqrt{Dt}}\right]}dt}

We can find \int_{0}^{T}{\text{erf}{\left[\frac{x}{2\sqrt{Dt}}\right]}dt}: First, letting u=t^{-\frac{1}{2}}, so \frac{du}{dt}=-\frac{1}{2t^{-\frac{3}{2}}\ }=-\frac{1}{2u^3}:

\int_{0}^{T}{\text{erf}{\left[\frac{x}{2\sqrt{Dt}}\right]}dt}=-2\int_{0}^{\frac{1}{\sqrt T}}{\frac{\text{erf}{\left[\frac{xu}{2\sqrt D}\right]}}{u^3}du\ }

Let v=\frac{u}{2\sqrt D} such that \frac{dv}{du}=\frac{1}{2\sqrt D} and u^{-3}=\frac{1}{8D^\frac{3}{2}v^3}: \int_{0}^{\frac{1}{\sqrt T}}{\frac{\text{erf}{\left[\frac{xu}{2\sqrt D}\right]}}{u^3}du\ }=\frac{1}{4D}\int_{0}^{\frac{1}{2\sqrt{DT}}}{\frac{\text{erf}{\left[xv\right]}}{v^3}dv\ }

Let p=\text{erf}{\left(xv\right)} and \frac{dq}{dv}=v^{-3} such that \frac{dp}{dv}=\frac{2xexp\left[-x^2v^2\right]}{\sqrt\pi} and q=-\frac{1}{2v^2}.

By parts:

\int_{0}^{\frac{1}{2\sqrt{DT}}}{\frac{\text{erf}{\left[xv\right]}}{v^3}dv=-\left.\frac{\text{erf}{\left[xv\right]}}{2v^2}\right|_0^{\frac{1}{2\sqrt{DT}}}+\int_{0}^{\frac{1}{2\sqrt{DT}}}\frac{x\exp{\left[-x^2v^2\right]}}{\sqrt\pi v^2}\ dv\ }\\ \int_{0}^{\frac{1}{2\sqrt{DT}}}\frac{xexp\left[-x^2v^2\right]}{\sqrt\pi v^2}\ dv=\frac{x}{\sqrt\pi}\int_{0}^{\frac{1}{2\sqrt{DT}}}\frac{\exp{\left[-x^2v^2\right]}\ }{v^2}\ dv

Let m=\exp{\left[-x^2v^2\right]} and \frac{dn}{dv}=v^{-2} such that \frac{dm}{dv}=-2x^2v\exp{\left[-x^2v^2\right]} and n=-v^{-1}.

By parts:

\int_{0}^{\frac{1}{2\sqrt{DT}}}\frac{\exp{\left[-x^2v^2\right]}\ }{v^2}\ dv=-\frac{\exp{\left[-x^2v^2\right]}}{v}-\int_{0}^{\frac{1}{2\sqrt{DT}}}{2x^2\exp{\left[-x^2v^2\right]}dv}

Let w = xv:

\int_{0}^{\frac{1}{2\sqrt{DT}}} 2x^2 \exp\left(-x^2 v^2\right) dv = \sqrt{\pi} x \int_{0}^{\frac{x}{2\sqrt{DT}}} \frac{2 e^{-w^2}}{\sqrt{\pi}} dw \\ = \sqrt{\pi} x \text{erf}\left( \frac{x}{2\sqrt{DT}} \right) \\ \Longrightarrow \frac{x}{\sqrt{\pi}} \int_{0}^{\frac{1}{2\sqrt{DT}}} \frac{\exp\left(-x^2 v^2\right)}{v^2} dv = -\frac{x \exp\left(-x^2 v^2\right)}{\sqrt{\pi} v} - x^2 \text{erf}\left( \frac{x}{2\sqrt{DT}} \right) \\ \Longrightarrow -\left. \frac{\text{erf}(xv)}{2v^2} \right|_0^{\frac{1}{2\sqrt{DT}}} + \int_{0}^{\frac{1}{2\sqrt{DT}}} \frac{x \exp\left(-x^2 v^2\right)}{\sqrt{\pi} v^2} dv = -\frac{\text{erf}\left( \frac{x}{2\sqrt{DT}} \right)}{2v^2} - \frac{x \exp\left(-x^2 v^2\right)}{\sqrt{\pi} v} - x^2 \text{erf}\left( \frac{x}{2\sqrt{DT}} \right) \\ \Longrightarrow \frac{1}{4D} \int_{0}^{\frac{1}{2\sqrt{DT}}} \frac{\text{erf}(xv)}{v^3} dv = -\left( \frac{\text{erf}\left( \frac{x}{2\sqrt{DT}} \right)}{8Dv^2} + \frac{x \exp\left(-x^2 v^2\right)}{4\sqrt{\pi} D v} + \frac{x^2 \text{erf}\left( \frac{x}{2\sqrt{DT}} \right)}{4D} \right) \\ = \int_{0}^{\frac{1}{\sqrt{T}}} \frac{\text{erf}\left( \frac{xu}{2\sqrt{D}} \right)}{u^3} du = -\frac{1}{2} \int_{0}^{T} \text{erf}\left( \frac{x}{2\sqrt{Dt}} \right) dt

Since v=\frac{u}{2\sqrt D} and u=\frac{1}{\sqrt t}:

\int_{0}^{T}{\text{erf}{\left[\frac{x}{2\sqrt{Dt}}\right]}dt}=\frac{\text{erf}{\left[\frac{xu}{2\sqrt D}\right]}}{u^2}+\frac{x^2\text{erf}{\left[\frac{xu}{2\sqrt D}\right]}}{2D}+\frac{xexp\left[-\frac{x^2u^2}{4D}\right]}{\sqrt{\pi D}u}\\ =\frac{x\sqrt T\exp{\left[-\frac{x^2}{4DT}\right]}}{\sqrt{\pi D}}+\left(T+\frac{x^2}{2D}\right)\text{erf}{\left[\frac{x}{2\sqrt{DT}}\right]}\\ \Longrightarrow\int_{0}^{T}{\varphi_s\left\{1-\text{erf}{\left[\frac{x}{2\sqrt{Dt}}\right]}\right\}dt}=\varphi_s\left(T-\left\{\frac{x\sqrt T\exp{\left[-\frac{x^2}{4DT}\right]}}{\sqrt{\pi D}}+\left(T+\frac{x^2}{2D}\right)\text{erf}{\left[\frac{x}{2\sqrt{DT}}\right]}\right\}\right)

Plugging this back into Eq. (5.4):

\therefore \varphi(x,t) = -N_b \mathcal{K} v_{\text{max}} \cdot \varphi_s\left( t - \left\{ \frac{x \sqrt{t} \exp\left( -\frac{x^2}{4Dt} \right)}{\sqrt{\pi D}} + \left( t + \frac{x^2}{2D} \right) \text{erf}\left( \frac{x}{2\sqrt{Dt}} \right) \right\} \right) \\ - DA \frac{\varphi_s}{\sqrt{\pi D}} \left\{ \frac{x \sqrt{\pi}}{\sqrt{D}} \text{erf}\left( \frac{x}{2\sqrt{Dt}} \right) - 2 \sqrt{t} \exp\left( -\frac{x^2}{4Dt} \right) \right\}

Numerically (where A is the cross-sectional area of the chamber ( 0.15m radius, so A\approx0.035m^2 ) (N_b was estimated in Eq. (2.9) as at minimum ~3.6\times{10}^{11}):

\varphi(x,t) = -3.6 \times 10^{11} \cdot 0.99977 \cdot 1.9 \times 10^{-19} \cdot 10.45 \left( t - \left\{ \left( 1.1 \times 10^4 \right) x \sqrt{t} \exp\left( -10^8 \frac{x^2}{t} \right) + \left( t + \left( 2 \times 10^8 \right) x^2 \right) \text{erf}\left( 10^4 \frac{x}{\sqrt{t}} \right) \right\} \right) \\ - 1.0 \times 10^{-5} \left\{ \left( 3.5 \times 10^4 \right) x \text{erf}\left( 10^4 \frac{x}{\sqrt{t}} \right) - 2 \sqrt{t} \exp\left( -10^8 \frac{x^2}{t} \right) \right\}

The validity of this approximation can be assessed graphically, where Desmos (the graphing calculator used) computes the integrals numerically and thus compare them. The orange curve is g\left(t\right)=\int_{0}^{t}f\left(x\right)dx where f\left(x\right)=1-\erf{\left[x\right]} and the red curve is h\left(t\right)=\int_{0}^{t}{\frac{f\left(x\right)^2}{K_m+f\left(x\right)}dx}:

Concentration of N2O with respect to time, given our estimated integral model, has been graphed in orange. In red, the exact integral has been graphed.

Given the proximity of the graphs, the approximation seems valid. On a large scale, the difference is more visible, yet still quite negligible, contributing to a difference of ~0.8%. Moreover, looking around the value of y=K_m, where we expect the two curves to behave most differently, we see instead that they remain roughly parallel, the difference between them seemingly constant

Figure 2, zoomed in on figure 1
Zooming in to the previous figure, showing the difference between the graphs as when t = 1.4 x 10-5. The graphs look parallel, and the error is of only 0.8%.

Graphing the L2 Norm as y=\int_0^t\of\begin|g(t)-h(t)|^2\ dx〗, the difference in values between the two integrals can be found for any x=\frac{x_0}{2\sqrt{Dt_0}}:

Figure 3, difference in values
The difference in values between the integrals as distance from the centre and time vary.

N2O concentration at the centre of the cell

We analyse the concentration of N2O at the centre of a cell, which we take to be x=1.2\times{10}^{-6}m. (13)

Taking the upper bound for v_{max}, our analytic solution reads.

\varphi\left(1.2\times{10}^{-6},t\right)=-3.6\times{10}^{11}\times0.99977\times1.9\times{10}^{-19}\times10.45\left(t-\left\{\left(1.3\times{10}^{-2}\right)\sqrt t\exp{\left[-\frac{1.4\times{10}^{-4}}{t}\right]}+\left(t+2.9\times{10}^{-4}\right)\erf{\left[\frac{0.012}{\sqrt t}\right]}\right\}\right)-1.0\times{10}^{-5}\left\{0.042\erf{\left(\frac{0.012}{\sqrt t}\right)}-2\sqrt t\exp{\left[-\frac{1.44\times{10}^{-4}}{t}\right]}\right\}

Our graphical solution looks as follows:

Figure 4, nitrous oxide concentrations at centre of cell, upper bound.
Variation in concentration of N2O at the centre of the cell with respect to time for the upper bound of Vmax.

This model thus predicts that with the given bacterial population of 3.6\times{10}^{11} - an estimated minimum which considers biosafety - a square-base box with side length 15cm and height 0.3m (chosen for practicality and convenience in a hospital environment), a batch of 0.0705mol of N2O will be broken down in /~13.1 minutes. A standardised full gas canister of height 250mm and diameter 203mm, would contain 60%\times\frac{2.5dm\times{(1.0dm)}^2\pi}{24mol/dm^3}\approx0.20mol of N2O.
To ensure the same maximised time for full N2O processing, we predict the box should have a height of \frac{0.2}{10.45}\times\frac{2}{{0.15}^2\pi}\approx0.85m

Taking the lower bound for v_{max}, our analytic solution reads. \varphi\left(1.2\times{10}^{-6},t\right)=-3.6\times{10}^{11}\times0.99977\times6.5\times{10}^{-20}\times10.45\left(t-\left{\left(1.3\times{10}^{-2}\right)\sqrt t\exp{\left[-\frac{1.4\times{10}^{-4}}{t}\right]}+\left(t+2.9\times{10}^{-4}\right)\erf{\left[\frac{0.012}{\sqrt t}\right]}\right}\right)-1.0\times{10}^{-5}\left{0.042\erf{\left(\frac{0.012}{\sqrt t}\right)}-2\sqrt t\exp{\left[-\frac{1.44\times{10}^{-4}}{t}\right]}\right} Our graphical solution looks as follows:

Nitrous oxide concentrations at centre of cell, lower bound
Variation in concentration of N2O at the centre of the cell with respect to time for the lower bound of Vmax.

This suggests instead a time of ~111 minutes for the same 0.0106mol of N2O - incredibly different. To account for this large range of uncertainty, an average was taken and will be quoted as the predicted average time taken to process one canister of N2O: 62 minutes.

This is of course a large draw-back to this model, which occurred because of a lack of data available to the translation rate of N2OR. The most useful step for a future team to take, experimentally, or otherwise, would be to determine the ribosome density (ribden) for N2OR in E. coli, which would greatly improve this model.

For Future Teams: How to Calculate Ribosome Density

This is one of our main modelling suggestions to future teams. The time taken is the second point at which [N2O] = 0: -N_b\mathcal{K}v_{max}\bullet\varphi_s\left(t-\left{\frac{x\sqrt t\exp{\left[-\frac{x^2}{4Dt}\right]}}{\sqrt{\pi D}}+\left(t+\frac{x^2}{2D}\right)\erf{\left[\frac{x}{2\sqrt{Dt}}\right]}\right}\right)-DA\frac{\varphi_s}{\sqrt{\pi D}}\left{\frac{x\sqrt\pi}{\sqrt D}\erf{\left(\frac{x}{2\sqrt{Dt}}\right)}-2\sqrt t\exp{\left[-\frac{x^2}{4Dt}\right]}\right}=0

Where: v_{max}=\eta\bullet\frac{\left[mRNA\right]\bullet r i b d e n\bullet r i b o c c}{K_m+\left[mRNA\right]\bullet r i b d e n\bullet r i b o c c}

\Longrightarrow\frac{ribden}{K_m+\left[mRNA\right]\bullet r i b d e n\bullet r i b o c c}=-\frac{DA\frac{\varphi_s}{\sqrt{\pi D}}\left\{\frac{x\sqrt\pi}{\sqrt D}\erf{\left(\frac{x}{2\sqrt{Dt}}\right)}-2\sqrt t\exp{\left[-\frac{x^2}{4Dt}\right]}\right\}}{N_b\mathcal{K}\varphi_s\eta\left[mRNA\right]\bullet r i b o c c\left(t-\left\{\frac{x\sqrt t\exp{\left[-\frac{x^2}{4Dt}\right]}}{\sqrt{\pi D}}+\left(t+\frac{x^2}{2D}\right)\erf{\left[\frac{x}{2\sqrt{Dt}}\right]}\right\}\right)} ribden=-\frac{\left[K_m\frac{DA\frac{\varphi_s}{\sqrt{\pi D}}\left\{\frac{x\sqrt\pi}{\sqrt D}\erf{\left(\frac{x}{2\sqrt{Dt}}\right)}-2\sqrt t\exp{\left[-\frac{x^2}{4Dt}\right]}\right\}}{N_b\mathcal{K}\varphi_s\eta\left[mRNA\right]\bullet r i b o c c\left(t-\left\{\frac{x\sqrt t\exp{\left[-\frac{x^2}{4Dt}\right]}}{\sqrt{\pi D}}+\left(t+\frac{x^2}{2D}\right)\erf{\left[\frac{x}{2\sqrt{Dt}}\right]}\right\}\right)}\right]}{1+\left[\frac{DA\frac{\varphi_s}{\sqrt{\pi D}}\left\{\frac{x\sqrt\pi}{\sqrt D}\erf{\left(\frac{x}{2\sqrt{Dt}}\right)}-2\sqrt t\exp{\left[-\frac{x^2}{4Dt}\right]}\right\}}{N_b\mathcal{K}\varphi_s\eta\left(t-\left\{\frac{x\sqrt t\exp{\left[-\frac{x^2}{4Dt}\right]}}{\sqrt{\pi D}}+\left(t+\frac{x^2}{2D}\right)\erf{\left[\frac{x}{2\sqrt{Dt}}\right]}\right\}\right)}\right]}

Experimentally, several points can be taken from some position x into the bacterial cell and the time t from the introduction of aqueous N2O into the system.

Important assumptions:

  1. This model assumed that all the bacteria are all in one layer covering the cross-sectional area of the chamber in which the breakdown takes place, and that the N2O enters from underneath and spreads evenly out to all of them.
  2. For this to be roughly equivalent to the current set-up, the N2O must be evenly mixed throughout the system such that each molecule can diffuse into a bacterium each at the same distance from it. This will ultimately be inaccurate, but to reduce this, mixing in the final cycle of design for the hardware has been optimised (see Design).
  3. This model assumes all N2O entering the system and dissolves simultaneously, though this would take a few seconds in real life. This effect was deemed negligible for this model, though we would invite future teams to investigate the effect this may have on the model.

References

(0) NHS Royal Devon University Healthcare, NHS Foundation Trust. Entonox. Leaflet number: 451. Version number: 2. p1. Available at: https://www.royaldevon.nhs.uk/media/giwp330t/entonox.pdf [accessed 5/9/24]

(1) Hendrickx, J.F.A., Cardinael, S., Carette, R., Lemmens, H. J.M., De Wolf, A.M. (2007). The ideal oxygen/nitrous oxide fresh gas flow sequence with the Anesthesia Delivery Unit machine. Journal of Clinical Anesthesia, 19(4), pp. 274-279. doi: https://doi.org/10.1016/j.jclinane.2007.01.003

(2) Sander, R. (2015). Compilation of Henry’s law constants (version 4.0) for water as solvent. Atmospheric Chemistry and Physics, 15, pp. 4399-4981. doi: https://doi.org/10.5194/acp-15-4399-2015

(3) Pinnick, E.R., Erramilli, S., Wang, F. (2010). The potential of mean force of nitrous oxide in a 1,2-dimyristoylphosphatidylcholine lipid bilayer. Chemical Physics Letters, 489(1-3), pp. 96-98. doi: https://doi.org/10.1016/j.cplett.2010.02.047

(4) Chen, J., Yue, K., Shen, L., Zheng, C., Zhu, Y., Han, K., Kai, L. (2023). Aquaporins and CO2 diffusion across biological membrane. Frontiers in Physiology, (14) doi: https://doi.org/10.3389%2Ffphys.2023.1205290

(5) Wang, S., Zhou, T., Trusler, M.J.P. (2023). Diffusion Coefficients of N2O and H2 in Water at Temperatures between 298.15 and 423.15 K with Pressures up to 30 MPa. Journal of Chemical Engineering Data, 68(6), pp. 1313-1319 doi: https://doi.org/10.1021/acs.jced.3c00085

(6) Simon, S.A., Gutknecht, J. (1980). Solubility of carbon dioxide in lipid bilayer membranes and organic solvents. Biochimica et Biophysica Acta, 596(3), pp. 352-8

doi: https://doi.org/10.1016/0005-2736(80)90122-4 (7) Cubaud, T., Sauzade, M., Sun, R. (2012). CO2 dissolution in water using long serpentine microchannels. Biomicrofluids, 6(2), pp. 22002-22009

doi: https://doi.org/10.1063%2F1.3693591 (8) Zhao, F., Li, R., Liu, Y., Chen, H. (2022). Perspectives on lecithin from egg yolk: Extraction, physicochemical properties, modification, and applications. Frontiers in Nutrition, 9. doi: https://doi.org/10.3389%2Ffnut.2022.1082671

(9) Rowlett, V.W., Mallampalli, V.K.P.S, Karlstaedt, A., Dowhan, W., Taegtmeyer, H., Margolin, W., Vitrac, H. (2017). Impact of Membrane Phospholipid Alterations in Escherichia coli on Cellular Function and Bacterial Stress Adaptation. Journal of Bacteriology, 199 (13), e00136-17 doi: https://doi.org/10.1128%2FJB.00849-16

(10) Moeller, M.N., Li, Q., Chinnaraj, M., Cheung, H.C., Lancaster Jr, J.R., Denicola, A. (2016). Solubility and diffusion of oxygen in phospholipid membranes. Biochimica et Biophysica Acta- Biomembranes, 1858(11), pp. 2923-2930 doi: https://doi.org/10.1016/j.bbamem.2016.09.003

(11) Wei, W., Mandin, C., Blanchard, O., Mercier, F., Pelletier, M. (2016). Temperature dependence of the particle/gas partition coefficient: An application to predict indoor gas-phase concentrations of semi-volatile organic compounds. Science of the Total Environment, 563–564, pp.506-512 doi: https://doi.org/ff10.1016/j.scitotenv.2016.04.106

(12) Park, M., Yoo, G., Bong, J-H., Joachim, J., Kang, M-J., Pyun, J-C. (2015). Isolation and characterization of the outer membrane of Escherichia coli with autodisplayed Z-domains. Biochimica et Biophysica Acta- Biomembranes, 1848(3), pp. 842-847 doi: https://doi.org/10.1016/j.bbamem.2014.12.011

(13) Yamada, H., Yamaguchi, M., Shimizu, K., Murayama, S.Y., Mitarai, S., Sasakawa, C., Chibana, H. (2017). Structome analysis of Escherichia coli cells by serial ultrathin sectioning reveals the precise cell profiles and the ribosome density. Microscopy, 66(4), pp. 283-294 doi: https://doi.org/10.1093/jmicro/dfx019

(14) The general solution was found, and the following derivations are framed around the working shown in notes from: Prof. Virkar, A. Lecture 4: Diffusion (Fick’s second law), 5034/6034: Chemical Engineering. John and Marcia Price College of Engineering

(15) Yan, Q., Fong, S.S. (2017). Study of in vitro transcriptional binding effects and noise using constitutive promoters combined with UP element sequences in Escherichia coli. Journal of Biological Engineering, 11(33) doi: https://doi.org/10.1186%2Fs13036-017-0075-2

(16) Chan, J., Davis, C., Jokic, I. (2006). Influences of growth temperature and preparation of competent cells on efficiency of chemically-induced transformation in Escherichia coli DH5α. Journal of Experimental Microbiology and Immunology, 9, pp. 92-96 Pdf accessed on [14/9/24] at: https://www.microbiology.ubc.ca/sites/default/files/roles/drupal_ungrad/JEMI/9/9-92.pdf

(17) Volkmer, B., Heinemann, M. (2011). Condition-Dependent Cell Volume and Concentration of Escherichia coli to Facilitate Data Conversion for Systems Biology Modeling. PLoS One, 6(7), e23126 doi: https://doi.org/10.1371%2Fjournal.pone.0023126

(18) Shepherd, N., Dennis, P., Bremer, H. (2001). Cytoplasmic RNA Polymerase in Escherichia coli. Journal of Bacteriology, 183(8), pp. 2527-2534 doi: https://doi.org/10.1128%2FJB.183.8.2527-2534.2001

(19) Bernstein, J.A., Khodursky, A.B., Lin, P-H., Cohen, S.N. (2002). Global analysis of mRNA decay and abundance in Escherichia coli at single-gene resolution using two-color fluorescent DNA microarrays. Proceedings of the National Academy of Sciences, 99(15), pp. 9697-9702 doi: https://doi.org/10.1073/pnas.112318199

(20) Brockmann, R., Beyer, A., Heinisch, J.J., Wilhem, T. (2007). Posttranscriptional Expression Regulation: What Determines Translation Rates? PLoS Computational Biology. Available at: https://doi.org/10.1371/journal.pcbi.0030057 [accessed 14/9/24]

(21) Beyer, A., Hollunder, J., Nasheuer, H.P., Wilhelm, T. (2004). Post-transcriptional Expression Regulation in the Yeast Saccharomyces cerevisiae on a Genomic Scale. Molecular and Cellular Proteomics, 3(11), pp. 1083-1092 doi: https://doi.org/10.1074/mcp.M400099-MCP200

(22) Matamouros, S., Gensch, T., Cerff, M., Sachs, C.C., Abdollahzadeh, I., Hendrik, J., Horst, L., Tenhaef, N., Tenhaef, J., Noack, S., Garf, M., Takors, R., Noeh, K., Bott, M, (2023). Growth-rate dependency of ribosome abundance and translation elongation rate in Corynebacterium glutamicum differs from that in Escherichia coli. Nature Communications, 14, 5611 doi: https://doi.org/10.1038/s41467-023-41176-y

(23) Dai, X., Zhu, M., Warren, M., Balakrishnan, R., Patsalo, V., Okano, H., Williamson, J.R., Fredrick, K., Wang, Y-P., Hwa, T. (2016). Reduction of translating ribosomes enables Escherichia coli to maintain elongation rates during slow growth. Nature Microbiology, 12(2), p. 16231 doi: https://doi.org/10.1038/nmicrobiol.2016.231

(24) Nguyen, H.L., Duviau, M.P., Laguerre, S., Nouaille, S., Cocaign-Bousquet, M., Girbal, L. (2022). Synergistic Regulation of Transcription and Translation in Escherichia coli Revealed by Codirectional Increases in mRNA Concentration and Translation Efficiency. Microbiology Spectrum, 10(1) doi: https://doi.org/10.1128%2Fspectrum.02041-21

(25) Carreira, C.C.S. (2017). Insights into the structure and reactivity of the catalytic site of nitrous oxide reductase. Ph. D. Thesis. Universidade Nova de Lisboa. Available at: https://run.unl.pt/bitstream/10362/27900/1/Carreira_2017.pdf [accessed 16/9/24]

(26) Liebecq, C. (1992). Biochemical Nomenclature and Related Documents. 2nd edition. Portland, Maine: Portland Press, p. 100 Available at: https://iubmb.qmul.ac.uk/kinetics/ek4t6.html

(27) Wilks, J.C., Slonczewski, J.L. (2007). pH of the Cytoplasm and Periplasm of Escherichia coli: Rapid Measurement by Green Fluorescent Protein Fluorimetry. Journal of Bacteriology, 189(15), pp. 5601-5607 doi: https://doi.org/10.1128%2FJB.00615-07

(28) Pan, Y., Ni, B-J., Yuan, Z. (2012). Effect of pH on N2O reduction and accumulation during denitrification by methanol utilizing denitrifiers. Water Research, 46(15), pp. 4832-4840 doi: https://doi.org/10.1016/j.watres.2012.06.003

(29) Lewis Jr., C.A., Wolfenden, R. (2008). Uroporphyrinogen decarboxylation as a benchmark for the catalytic proficiency of enzymes. Proceedings of the National Academy of Sciences, 105(45), pp. 17328-17333

(30) Yi, L., Sheng, R., Wei, W., Zhu, B., Zhang, W. (2022). Differential contributions of electron donors to denitrification in the flooding-drying process of a paddy soil. Applied Soil Ecology, 177, 104527 doi: https://doi.org/10.1016/j.apsoil.2022.104527

(31) Pinochet-Barros, A., Helmann, J.D. (2018). Redox Sensing by Fe2+ in Bacterial Fur Family Metalloregulators. Antioxidants and Redox Signaling, 29(18), pp. 1858-1871 doi: https://doi.org/10.1089%2Fars.2017.7359

(32) Martin, J.E., Waters, L.S., Storz, G., Imlay, J.A. (2015). The Escherichia coli Small Protein MntS and Exporter MntP Optimize the Intracellular Concentration of Manganese. PLoS Genetics, 11(3): e1004977 doi: https://doi.org/10.1371%2Fjournal.pgen.1004977

(33) Zhang, L., Wust, A., Prasser, B., Muller, C., Einsle, O. (2019). Functional assembly of nitrous oxide reductase provides insights into copper site maturation. Proceedings of the National Academy of Sciences, 116(26), pp. 12822-12827 doi: https://doi.org/10.1073%2Fpnas.1903819116

(34) Barnaba, C., Taylor, E., Brozik, J.A. (2017). Dissociation Constants of Cytochrome P450 2C9 / Cytochrome P450 Reductase Complexes in a Lipid Bilayer Membrane Depend on NADPH: A Single Protein Tracking Study. Journal of the American Chemical Society, 139(49), pp. 17923-17934 doi: https://doi.org/10.1021/jacs.7b08750

(35) Rathnayaka, S.C., Mankad, N.P. (2021). Coordination chemistry of the CuZ site in nitrous oxide reductase and its synthetic mimics. Coordination Chemistry Reviews, 429, 213718 doi: https://doi.org/10.1016/j.ccr.2020.213718

(36) Ebert, A., Hannesschlaeger, C., Goss, K-U., Pohl, P. (2018). Passive Permeability of Planar Lipid Bilayers to Organic Anions. Biophysical Journal, 115(10), pp.1931-1941 doi: https://doi.org./10.1016/j.bpj.2018.09.025

Section 2: ComC - ComR evaluation

Evaluation

It is still unclear how copper enters the bacterial cytoplasm. In E. coli, permeability to copper is reduced by the ComC outer membrane protein through an unknown mechanism. ComC expression decreases under low copper conditions to allow copper uptake. Therefore, unfortunately there is no way yet to mathematically find a relationship between the specific effect of ComC on permeability to copper ions.

However, there have been studies conducted to experimentally verify the effect of ComC on copper uptake:(1)

Graph showing relative luminescense against CuSO4 concentration

Luminescence was recorded quantitatively (using a lux based biosensor) in liquid media for a 𝚫ComR (Reduced ComR) strain (open circle) and the wild-type (closed circle) at different extracellular CuSO4 concentrations. Luminescence was normalized to the number of cells. The reduced luminescence of the 𝚫ComR strain suggests that this mutant contained lower cytoplasmic copper levels compared to the wild-type when exposed to copper.

Similarly, below is a graph (1) demonstrating the expression of the ComR (open circle) and the ComC (closed circle) genes in wild-type E. coli, determined by qPCR at a range of CuSO4 concentrations.

Graph showing relative mRNA levels against CuSO4 concentration

Expression of the ComR (open circle) and the ComC (closed circle) gene in wild-type E. coli, determined by qPCR at a range of CuSO4 concentrations.

The expression levels of ComR and ComC in response to copper were assessed by qPCR (to measure gene expression). Expression of ComR was not significantly stimulated by copper. In contrast, ComC was induced approximately 30-fold by 3 mM copper. The expression of ComC was also determined in the ComR knockout strain and was found to be increased by 270-fold on average, showing the repression of ComC by ComR in the wild-type.

Two graphs showing relative cus and copA mRNA levels with and without copper

To assess periplasmic copper levels, the study made use of the CusRS regulatory system which naturally senses periplasmic copper and regulates the cusCFBA operon in response to it. When the 𝚫ComC mutant was exposed to copper, expression of the cusCFBA operon was induced over 50% more strongly than in the wild-type under the same conditions, as assessed by qPCR. This points to increased copper leakage across the outer membrane in the 𝚫ComC mutant. As expected, this in turn also leads to increased cytoplasmic copper, as evident by the increased expression of copA, which is under the control of the cytoplasmic copper-sensing CueR regulator.(2)

Therefore, we can see that ComR represses ComC, a gene encoding an outer membrane protein that reduces copper permeability in E. coli. The absence of ComC results in higher copper levels in both the periplasm and cytoplasm, increasing copper sensitivity. These results gave us confidence in the premise behind our project.

Furthermore, gauging the intricacies of the interaction between ComR and ComC proves to be challenging due to there being a gap in experimental data for the precise nature of their relationship. In order to model the interaction between the two proteins, it is possible to deploy molecular dynamics simulations as a method of computational protein-protein analysis as it would allow for a predicted result for the likelihood and detail of the interaction.

The software PEPPI (Pipeline for the Extraction of Predicted Protein-protein Interactions) by The Zhang Lab Group (3) is one such example of a resource that was implemented when modeling ComR and ComC. However, due to the novelty of both proteins’ properties in regards to each other, it is evident that this computational approach requires further investigation for a succinct and meticulous numerical explanation.

Therefore, for the purposes of the design, the assumption made for ComR and ComC are that they follow a 1:1 relationship when interacting i.e. upregulation of one unit of ComR will result in a downregulation of one unit of ComC.

References

(1) Mermod, M., Magnani, D., Solioz, M. and Stoyanov, J.V. (2011). The copper-inducible ComR (YcfQ) repressor regulates expression of ComC (YcfR), which affects copper permeability of the outer membrane of Escherichia coli. BioMetals, 25(1), pp.33–43. doi:https://doi.org/10.1007/s10534-011-9510-x.

(2) Porcheron, G., Garénaux, A., Proulx, J., Sabri, M. and Dozois, C.M. (2013). Iron, copper, zinc, and manganese transport and regulation in pathogenic Enterobacteria: correlations between strains, site of infection and the relative importance of the different metal transport systems for virulence. Frontiers in Cellular and Infection Microbiology, 3. doi:https://doi.org/10.3389/fcimb.2013.00090. Microbiol. 2013 Dec 5;3:90. doi: 10.3389/fcimb.2013.00090. PMID: 24367764; PMCID: PMC3852070.

(3) https://zhanggroup.org/PEPPI/

Section 3: Modelling implementation costs

Volume in canisters / year

Let us begin by defining the number of maternity deliveries in the NHS as N_{births}. N_{births} = 547, 224 in 2022-2023 (1). Let us now define the proportion of women choosing to use Entonox when giving birth as P_{use}, which an NHS 2022 Survey reveals stands at 0.76 (2).

As such, N_{use} = N_{births} \times P_{use}

We now figure out the total volume exhaled by women into canisters, to figure out the total volume of gas that will be treated per year. Women breathe in Entonox intermittently for each contraction, meaning they will take a few breaths of Entonox, followed by breaths of regular air, as it starts to relieve pain approximately 20 seconds after inhalation (3). Between contractions, women do not breathe in Entonox, and thus, our design does not collect air exhaled between contractions.

Total volume exhaled in a contraction may be written as:

V_{contraction} = V_{tidal} \times B_{contraction}

Where we take the tidal volume of a breath during labor to be 750 ml or 0.00075m^3 (4), and the number of breaths per contraction B_{contraction}, to be 34.2 (4).

We can figure out the total number of contractions, N_{contraction}, by using the formula:

N_{\text{contraction}} = \frac{T_{\text{Stage I}}}{T_{\text{contraction}} + T_{\text{rest}}}

Where T_{Stage I} is the average labor time spent in Stage I labor, with regular contractions, in hospitals. Women are typically admitted into hospitals when their cervix is dilated from at 4cm (5), and progress into Stage II labor, when their cervix is dilated at 10cm. The American College of Obstetricians and Gynecologists (ACOG) have noted that the average dilatation rate is 1cmh-1, meaning that women spend an average of 6 hours in Stage I labour in NHS hospitals. As such, we take T_{Stage I} to be 21600s. We take T_{contraction}, the time spent in each contraction, to be 46s (4), and T_{rest}, the average time spent in between contractions, to be 3 minutes, or 180s (6).

Combining these equations, we can find the total volume, V_{tot}, to be processed in a year to be

V_{\text{tot}} = N_{\text{use}} \times V_{\text{contraction}} \times N_{\text{contraction}} \\ V_{\text{tot}} = N_{\text{births}} \times P_{\text{use}} \times V_{\text{tidal}} \times B_{\text{contraction}} \times \frac{T_{\text{Stage I}}}{T_{\text{contraction}} + T_{\text{rest}}}

The assumptions we make are as follows:

  1. We recognize that each labor process is vastly different—and thus, the volume exhaled by each woman will be vastly different. We hope that these values average out to the ones we have used, given the large number of births per year. In particular, we assume that the capacity of lungs, number of breaths exhaled into the mouthpiece into the canister, length of labor, rest-time between contractions, and contractions themselves are similar.
  2. (ii) In particular, we assume that the capacity of lungs, number of breaths exhaled into the mouthpiece into the canister, length of labor, rest-time between contractions, and contractions themselves are similar.

Cost of machines running

In a year, N_{\text{canisters}} will be produced, where:

N_{\text{canisters}} = \frac{V_{\text{tot}}}{V_{\text{canister}}}

We can find V_{\text{canister}} by considering a cylindrical, standardized, full gas canister which would have a radius of 101.5 mm or 0.1015 m and a height of 250 mm or 0.25 m.

V_{\text{canister}} = \pi \left( R_{\text{canister}} \right)^2 H_{\text{canister}}

As such, N_{\text{canisters}} = \frac{V_{\text{tot}}}{\pi \left( R_{\text{canister}} \right)^2 H_{\text{canister}}}

From our N_2O reduction modeling, we saw that each canister would take 64 minutes, or 1.067 hours, to be processed. Total time T to process all canisters would be:

T = T_{\text{canister}} \times N_{\text{canisters}}

Let us now take a conservative estimate that our machines must undergo maintenance for 4 hours per day, and thus run only 20 hours per day. The total number of machines, N_{\text{machines}}, would be:

N_{\text{machines}} = \frac{T}{365 \times 20}

The total cost of machines would need to factor in the cost per unit, which we calculated to be approximately £160. The total cost of machines can be written as:

C_{\text{total}} = C_{\text{machine}} \times N_{\text{machines}}

To combine all these separate parts into a large equation:

C_{\text{total}} = \frac{C_{\text{machine}} \cdot T_{\text{canister}} \cdot N_{\text{births}} \cdot P_{\text{use}} \cdot V_{\text{tidal}} \cdot B_{\text{contraction}} \cdot T_{\text{Stage I}}}{365 \cdot 20 \cdot \pi \cdot \left( R_{\text{canister}} \right)^2 \cdot H_{\text{canister}} \cdot \left( T_{\text{contraction}} + T_{\text{rest}} \right)}

Substituting all the relevant values:

C_{\text{total}} = \frac{160 \cdot 1.067 \cdot 547224 \cdot 0.76 \cdot 0.00075 \cdot 34.2 \cdot 21600}{365 \cdot 20 \cdot \pi \cdot \left( 0.1015 \right)^2 \cdot \left( 0.25 \right) \cdot \left( 46 + 180 \right)}

The total cost C_{\text{total}} stands at roughly £2,947,000. In short, it would cost roughly £2,947,000 to build enough machines to treat all Entonox emissions from maternity units in the first year. Once these machines are built, we would hope they would not need to be replaced later. These costs exclude the energy costs from running the machines, transportation costs of canisters to offshore factories containing the machines, and staff employed to run the machines, which would further increase the cost.

As a reference point, the total NHS funding in 2022-23 stood at £153 billion.

References

(1) NHS Maternity Statistics, England, 2022-23. Accessed 30/09/2024. Link: https://digital.nhs.uk/data-and-information/publications/statistical/nhs-maternity-statistics/2022-23

(2) Drew, Rebecca. Women in labor denied Entonox because of poor ventilation in delivery suites. Fieldfisher. 23/03/2023. Link: https://www.fieldfisher.com/en/injury-claims/insights/women-in-labour-denied-entonox-because-of-poor-ventilation

(3) A. Young, M.Ismail, AG Papatsoris et. Al. Entonox inhalation to reduce pain in common diagnostic and therapeutic outpatient urological procedures: a review of the evidence. Ann R Coll Surg Engl. 2012 Jan; 94(1): 8-11. DOI: 10.1308/003588412X13171221499702

(4) J. Selwyn Crawford, M.E. Tunstall. Notes on respiratory performance during labour. Brit. J. Anaesth (1968), 40, 612. PDF: https://www.bjanaesthesia.org.uk/article/S0007-0912(17)52195-9/pdf

(5) The stages of labour and birth. NHS. Accessed 30/09/2024. Link: https://www.nhs.uk/pregnancy/labour-and-birth/what-happens/the-stages-of-labour-and-birth/

(6) Stages of Labor. Cleveland Clinic. Accessed 30/09/2024. Link: https://my.clevelandclinic.org/health/symptoms/22640-stages-of-labor

(7) Our funding. NHS England. Accessed 30/09/2024. Link: https://www.england.nhs.uk/publications/business-plan/our-2022-23-business-plan/our-funding/