Quorum Sensing

For the final product, BeeT, we intend to use toxins produced by Bacillus thuringiensis, also called Cry toxins. However, these toxins are also harmful to our chassis, resulting in a reduction of toxin production when it is needed to kill the mites. To counteract this effect we envision the use of quorum sensing that activates Cry toxin production only when there is a large quantity of BeeT present. Ideally, BeeT is able to produce Cry toxin over a long period to improve its effectiveness against mites. However, if an entire population of BeeT is synchronized, we hypothesize that only a single burst of Cry toxin will take place before both the BeeT and mites are killed. Thus, this system will not be maximally effective over long time periods. To accomplish this we need multiple subpopulations of BeeT: Some producing BeeT while others are recuperating.To help understand this complex system we use dynamic modelling.

Genetic circuit

The quorum sensing system consists of two proteins LuxI and LuxR. LuxR is constitutively expressed, together with Acylated Homoserine Lactone (AHL) it forms a complex which activate the transcription of GFP ². AHL is a molecule that can diffuse freely through the cell membrane, and in this way diffuses from cell-to-cell. In a natural system the complex LuxR-AHL controls transcription of LuxI, this is a positive feedback loop that increases the amount of AHL in the system ³. When there is more AHL in the system, AHL is more likely to bind to the LuxR protein. With the use of quorum sensing we tried to simulate these different populations. In figure 1 you can see a schematic representation of the system.

What is Quorum Sensing?

Quorum sensing is a cell-cell communication system. The detection of chemical molecules allows the bacteria to distinguish between low and high cell densities. In this way the bacteria control gene expression in response to changes in cell number ^{1 4}. This process is achieved through the production and release of an in our case AHL autoinducer An autoinducer is a molecule that can diffuse through the cell membrane. . The AHL can travel from one cell to the other. There are many different types of autoinducers in quorum sensing systems. When sensed, the autoinducer can trigger other cells to produce more autoinducers.

According to the team from Davidson College and Missouri Western State University 2011 ⁵,negative feedback is present when LuxR protein is present and AHL-3OC6 is absent. They discovered that the BBa_K199052 part promotes "backwards transcription". In our project the importance of negative feedback in the quorum sensing system, to create different cell populations, has been investigated. As you can see in Figure 2 different populations can occur when the production rate of luxR and the complex forming are changed. However, when we remove the feedback in the system and change the same parameters, we get similar responses of the system. Shown in Figure 3.

This shows that the feedback has a minimal influence on the system to create different cell populations. Also when there is no negative feedback in the system parameter sets, we found more parameter sets that give a high amount of GFP. See Figure 4 and 5.

The figures are simulated with the GFP response obtained from the quorum sensing model, using the parameter sets that produce the best response. In the presence and absence of negative feedback in, the system different populations can occur by assuming the production rate of luxR and the production rate of the complex are different in each cell.
This system would work if we could control the LuxR production rate, however it could not be envisioned. Therefore, we designed another system that could help us create subpopulations downstream of the quorum sensing mechanism. The two systems were tested in the lab parallel with the modeling.

Why include a subpopulation system?

The subpopulation system consist of two genes; the first encodes for the protein that inhibits the systems expression, the other encodes for the corresponding activation protein of the system. The 434-cl-LVA inactivates the λ-cl directly, or prevents the translation of the λ-cl protein. 434-cl-LVA has a higher turnover rate than λ-cl. This subpopulation system is based on the system Bokinsky uses ⁶.

As you can see in Figure 6 glucose has a suppressing function on the system. Arabinose has a activating function on the system. We used the model to predict what will happen when we add glucose in different amounts. The assumptions we made for the subpopulation system are the following; 434 degrade faster than λ-cl, cell division happens when the cell size has doubled, and two cell populations exist with one growing differently than to another. We tested a number of parameter sets and picked the parameter set that gave the biggest RFP response based on the strength of glucose. Within the heat map you can see in which ratios the initial amounts of λ and 434 in the system are needed to get high RFP production. Where glucose and arabinose have a fixed number.

If λ and 434 are present in similar amounts we have a subpopulation system. But if they are present in different amounts than there will be no subpopulations. This can be expected when you look at the subpopulation system. The system is inhibited by 434, which represses the RFP production, and λ activates the RFP production. In figure 4 you can see balance between the 434 and λ amounts that are present for the output of RFP. This means that the initial conditions do not have much influence on the 434 and λ. With this particular parameter set, we can conclude that the translation rates are more important for the RFP response than the initial conditions. For other parameter sets the initial conditions are important. To understand how changing the RBS impacts RFP responses, we simulated the system for different combinations of transcription rates. From these simulations we found:

To confirm that our model represents Biological behaviour we analysed the RBS library. We can conclude out of these results that the RFP quantitatively corresponds to the Figure 7

Combined system

Based on the modeling we hypothesized the following: When there are more cells are present in the system, more AHL-LuxR complex is formed. The complex inhibits the subpopulation promoter. When the promoter is inhibited production of 434 will be suppressed and production of λ cl will be activated. At a certain time point λ-cl takes control over the system, because 434 has a higher turnover rate than λ-cl. In this case more λ-cl results in more RFP.

In a later stage, to get the desired response, the quorum sensing system and the subpopulation system were combined. As shown in Figure 8, you can see how we think to generate different cell populations. With this extended part we try to generate a system which predicts the behaviour of the subpopulations. In figure 9 you can see the increasing RFP over time after many cell divisions, indicating an increasing cell population.

Methods

During the research Matlab version R2016a has been used.
Because there was no data from the wet lab we assumed that all parameters exist within a biologically reasonable range numbers between 0 and 1. For the analyses we used 100.000 parameter sets, 10 cells, random initial conditions and 60 timesteps. Tuneable parameters are used, each parameter set can produce dramatically different population dynamics. To determine which of these parameters produce the best system response we used Latin Hypercube Latin hypercube is a statistical method to get random numbers from a box of x by n numbers. For example, if x = 4, where x is the number of divisions within the parameter value range, and n = 2, where n is number of parameters, you will obtain a box with 4 square times 2 square, giving you 24 random numbers. Within each division of the parameter space a single random number is chosen. sampling.

With the parameter numbers of the Latin Hypercube sampling we made parameter sets that are obtained from a lognormal distribution a lognormal distribution assigns probabilities to all positive values. However the distributions is skewed to favour smaller values. Biological reaction rates have been found to follow this distribution 8 with a parameter estimation based on Raue et al ⁷. Below you could see an example;

With these confidence intervals Equation for confidence interval used the best parameter sets could be chosen.

In conclusion

Light Kill Switch

To prevent BeeT from escaping into the environment around the hive and spreading we built a light kill switch. This system consists of a light switch that triggers when blue light hits the organisms. This light switch is coupled to a toxin/anti-toxin system, when the light switch is triggered it inhibits the production of the anti-toxin. With the anti-toxin production inhibited, the constitutive production of toxin kills BeeT. This system was also modeled using dynamic modelling, this was to ensure that the system only kills BeeT in the presence of blue light and it survives when no blue light is present.

Light on, BeeT off

One important aspect of our project is the 1. The cell death is inflicted by the mazEF toxin-antitoxin (TA) system. The latter functions by mazE (antitoxin) preventing the toxin (mazF) from cleaving mRNAs which ultimately kills the cell.
With the underlying biology of the system in mind, we then used synthetic bioengineering principles to mathematically describe our design and explore the dynamics of the system.

Our results show that we have a working mathematical model which describes the data of Ohlendorf et al.¹. With this model we now know where to look for the relative kinetic rates in the parameter space to give suggestions to the wet-lab on how to build a functioning optogenetic kill switch. Our differential equations of pDusk and pDawn can also be used by future iGEM teams to mathematically describe and build new optogenetic switches. The Matlab code is made available here.

Light on - Modelling The Optogenetic Switch

For the first part of our analysis, we optimized the mathematical models of pDusk and pDawn using published data¹. For now, we ignored the toxin-antitoxin system. This gave us the option to construct our scoring function to evaluate both systems simultaneously with a weighted means approach.
With the construction of the optogenetic switch part of our overall kill switch system, we are confident to represent the behaviour of the microorganism in the lab mathematically. Suggestions from the model can then be used as a reference point for building the toxin-antitoxin system into the kill switch model.

Goals

• Parameter Estimation: search the parameter space to identify parameter sets which satisfy the scoring, i.e. describe the behaviour of the Figure 2a in Ohlendorf et al.¹.
• Model Code: provide a mathematical model which can be used by future iGEM teams to built optogenetic switches.

System Design

We first studied the pDusk system in more detail to derive simplified equations which can describe the system behaviour. The genetic switches are based on the interactions between YF1, the light sensing protein, and FixJ, the response regulator as described by Möglich et al.². In darkness, YF1 phosphorylates FixJ which in turn, activates the promoter pFixK₂. This activation leads to the expression of desired genes. Under light conditions, YF1 enters an excited light-state, which inhibits the phosphorylation of FixJ. In addition, it prevents the activation of the pFixK₂ promoter and thus represses the desired genes (see Figure_1a).

The pDawn system functions in the same way with the addition of cI_λ being expressed depending on the activity of the pFixK₂ promoter. The protein form of cI_λ then inhibits the expression of the RFP reporter gene (see Figure_1b).

Assumptions

To derive mathematical equations describing the above systems, we needed to make some simplifying assumptions.

• we neglect the effects of dilution and growth (i.e. a cell’s volume is constant over time).
• we assumed, that the components in our system are not affected by and do not affect other cellular mechanisms.
• as there is no additional information gained from modeling transcription and translation for the dark state of YF1(y_DD) and the inactive form of FixJ (j_i) we lumped those variables and parameters. In addition, less variables and parameters translate into a less computationally expensive system.
• we also assumed the production, binding and degradation of mRNAs and proteins to be linear which resulted in a first order differential equations system

To model the light influence on the system, we assumed light activation to occur similar to the model on phytochrome B dimerization described in Klose et al.³ ($N\cdot k_2$). In addition, we consider in our system the different stages of YF1 (y_DD, y_DL/LD, and y_LL) as described by Möglich et al.². Translated to our system, this means that with increasing light intensities (N), the pool of activated YF1 (y_DL/LD and y_LL) would deplete the dark state of YF1 (y_DD), which is responsible for the phosphorylation and thus activation from inactive FixJ (j_i) to active FixJ (j_a) (see Figure_3). As all stages of YF1 form a dimer, we assumed dimerization to occur very fast and be the dominant and relevant form for our system. Lastly, we assumed a quasi steady-state approximation (Michaelis-Menten kinetics) for the effect of activated FixJ (j_a) on the expression downstream of the RFP (RFP_m) reporter gene.

Results

With our mathematical model we can describe the behaviour of both systems as tested in the lab by Ohlendorf et al.¹ (see Figure_2). We can thus say that we built a simplified mathematical representation of the pDusk and pDawn system which can describe the system behaviour.

1. In both pDusk and pDawn the black dots ± standard deviation represent the fluorescence over OD600 of the reporter gene RFP as tested in the lab by Ohlendorf et al.¹ (see their Figure 2a). The blue circles represent the computed RFP responses per light intensity normalized to the highest light intensity in each graph. The red line is given to show the trend of the behaviour of the system.

Preliminary results from the lab show, that the RFP_p response after 13 [h] inoculation time does not work in pDusk. The experiment was conducted in darkness and under 1554.7 [$\mu mol\cdot m^{-2}\cdot h^{-1}$] light intensity (see Table 2). The control showed a fluorescence of 47.9 * 10³ [a.u.], whereas the blank gave a response of 61.2 * 10³ [a.u.].

System	Dark	Light
pDusk	$15.4\cdot {10}^3 [a.u.]$	$15.2\cdot {10}^3 [a.u.]$
pDawn	11.0 ⋅ 10 - 3   [ a . u . ]	61.4 ⋅ 10 - 3   [ a . u . ]

If we simulate these situations in the model, we see the response shown in Figure 3. Looking at the time point at 13 [h] we can see a clear decrease from the RFP response from more than 0.2 [$\mu \mathrm{mol}$] to less than 0.05 [$\mu \mathrm{mol}$]. As we could expected a higher change in RFP within the pDusk system, it can be concluded, that either the dark conditions were not fully dark, or a more in depth parameter analysis would need to be conducted to see potential causes of the discrepancies.

1.

We also made the code for both systems available for future iGEM teams to build new optogenetic tools: Model.

Methods

After studying the pDusk and pDawn system, we constructed a parameter estimation procedure by picking random parameter sets based on the latin hypercube Latin hypercube is a statistical method to get random numbers from a box of x by n numbers. For example, if x = 4, where x is the number of divisions within the parameter value range, and n = 2, where n is number of parameters, you will obtain a box with 4 square times 2 square, giving you 24 random numbers. Within each division of the parameter space a single random number is chosen. sampling (LHS) method in Matlab 2015b (used for the whole study). We simulated 10 000 parameter sets and assumed the initial conditions of only y_DD and j_a being in the on state, i.e. we assigned the value one. All other initial stages were kept at 0, as there was no reason to assume otherwise. In addition we needed to transform those parameters to a more biologically feasible distribution as LHS generates uniformly distributed parameters. For those parameters where there was no information available, we assumed a widely biologically applicable lognormal distribution⁷. For the two parameters where we had data, we transformed those according to their known distribution (see Table 1 k₂ and k₃. We then loaded the RFP response data ± standard deviation from Ohlendorf et al.¹ for both pDusk and pDawn and transformed the light intensities from $[\mu W\cdot c{m}^{-2}]$ to $[\mu \mathrm{mol}\cdot m^{-2}\cdot h^{-1}]$.

We then simulated both systems according to our ordinary differential equations (ODEs) with the same parameter (sub)-sets and stored the response at 17 h for later evaluation in the scoring function. We chose to simulate the system until 17 h, as this is the experimental length Ohlendorf et al.¹ used in their study. As a next step we normalized the computed RFP_p response with the highest light intensity data point of Figure 2a from Ohlendorf et al.¹.

To evaluate how well a certain parameter set describes the response of the Ohlendorf et al.¹ system, we used the sum of squared residuals to score each parameter set as described in Raue et al.⁴:

${\chi }^2\left(\theta \right)=\sum _{k=1}^m\sum _{l=1}^d{\left(\frac{y_{kl}^d-y_k(\theta ,t_l)}{{\sigma }_{kl}^d}\right)}^2$

Where $y_{kl}^d$ represents $d$ data points for each observable $k$ at time points $t_l$. ${\sigma }_{kl}^d$ are the corresponding measurement errors and yk(θ,tl) the $k$^th observable as predicted by parameters $\theta$ for time point $t_l$.

To combine a corresponding score from the same parameter set for pDusk and pDawn we used the weigthed means approach accordingly:

$\frac{{\zeta }_u+{\zeta }_a}{2}$

Where ${\zeta }_u$ represents the score of a parameter set for pDusk and ${\zeta }_a$ the corresponding score of the same parameter set for pDawn.

To put the equations and parameters into context, we created a more detailed graphical representation of the pDusk and pDawn system in Figure_3 and provide the ordinary differential equations for both systems below.

Ordinary Differential Equations for pDusk

$\frac{{dy}_{DD}}{dt}=k_1+2\cdot k_2\cdot y_{DL,LD}-2\cdot (N\cdot k_3)\cdot y_{DD} -{\beta }_1\cdot y_{DD}$

$\frac{{dy}_{DL,LD}}{dt}=2\cdot (N\cdot k_3)\cdot y_{DD}+2\cdot k_2\cdot y_{LL}-2\cdot k_2\cdot y_{DL,LD}-2\cdot (N\cdot k_3)\cdot y_{DL,LD}-{\beta }_2\cdot y_{DL,LD}$

$\frac{{dy}_{LL}}{dt}=2\cdot (N\cdot k_3)\cdot y_{DL,LD}-2\cdot k_2\cdot y_{LL}-{\beta }_3\cdot y_{LL}$

$\frac{d{j}_i}{dt}=k_4+k_5\cdot j_a-{\beta }_4\cdot j_i$

$\frac{d{j}_a}{dt}=k_6\cdot y_{DD}\cdot j_i-{\beta }_5\cdot j_a$

$\frac{dRF{P}_m}{dt}=\frac{V_{max}\cdot j_a}{K_m+j_a}-{\beta }_6\cdot RF{P}_m$

$\frac{dRF{P}_p}{dt}=k_7\cdot RF{P}_m-{\beta }_7\cdot RF{P}_p$

Ordinary Differential Equations for pDawn

For pDawn we used the same equations as for pDusk from y_DD to j_a and added the following differential equations.

$\frac{dc{I}_m}{dt}=\frac{V_{max}\cdot j_a}{K_m+j_a}-{\beta }_8\cdot c{I}_m$

$\frac{dc{I}_p}{dt}=k_8\cdot c{I}_m-{\beta }_9\cdot c{I}_p$

$\frac{dRF{P}_m}{dt}=k_9\cdot \frac{1}{1+{\left(\frac{c{I}_p}{K_d}\right)}^2} -{\beta }_6\cdot RF{P}_m$

$\frac{dRF{P}_p}{dt}=k_7\cdot RF{P}_m-{\beta }_7\cdot RF{P}_p$

Here we show the best parameter set obtained from the parameter estimation (see Table 1). This parameter set was used as the basis for further analysis implementing the mazEF toxin-antitoxin system. To better put these numbers in context we provide a detailed graphical system representation of pDusk and pDawn in Figure_3.

Parameter	Value	Description
k1	$2.6921\frac{\mu \mathrm{mol}}{h}$	production rate of yDD
k2	$0.0008\frac{1}{h}$	relaxation rate of yDL,LD and yLL. We assumed a search space for τ of 5900 ± 25 s based on data from Möglich et al.2. Which was mathematically transformed to $k_2=\frac{\log \left(2\right)}{\tau }\cdot 3600\frac{s}{h}$.
k3	$0.4219\frac{m^2}{\mathrm{\mu mol}}$	conversion cross-section σ of light intensity activated production rate of yDL,LD and yLL. The search space for this parameter was defined as 1000 ± 25034.
β1	$0.3049\frac{1}{h}$	degradation rate of yDD
β2	$0.8406\frac{1}{h}$	degradation rate of yDL,LD
β3	$0.1477\frac{1}{h}$	degradation rate of yLL
k4	$0.2040\frac{\mathrm{\mu mol}}{h}$	production rate of ji
k5	$2.1623\frac{1}{h}$	de-phosphorylation rate of ja
β4	$0.5205\frac{1}{h}$	degradation rate of ji
k6	$2.0838\frac{1}{\mathrm{h\; \times \; \mu mol}}$	production rate of ja depending on the concentration of yDD and ji
β5	$0.6615\frac{1}{h}$	degradation rate of ja
Vmax	$2.9063\frac{\mathrm{\mu mol}}{h}$	Vmax of production rate of RFPm based on ja
KM	$0.7130\mu \mathrm{mol}$	KM of production rate of RFPm based on ja
β6	$2.0224\frac{1}{h}$	degradation rate of RFPm
k7	$0.0460\frac{1}{h}$	translation rate from RFPm to RFPp
β7	$0.2903\frac{1}{h}$	degradation rate of RFPp
β8	$1.1579\frac{1}{h}$	degradation rate of lambda phage inhibitor RNAm
k8	$3.8073\frac{1}{h}$	production rate of cIp depending on cIm
β9	$0.6563\frac{1}{h}$	degradation rate of cIp
k9	$0.9920\frac{\mathrm{\mu mol}}{h}$	maximal production rate of RFPm
KD	$0.1384\mathrm{\mu mol}$	dissociation constant of cIp at RFPm promoter. The Hill coefficient was chosen to be 2 as the cIp regulated promoter BBa_R0051 has 2 binding sites for cIp.

Conclusion

With the abovementioned results, we conclude to have built a mathematical model of pDusk and pDawn which describes the system behaviour. We know now where in the parameter space we would need to zoom in to define relative kinetic rates in collaboration with experimentalists. Initial results from the lab, did not show a response in the pDusk system, however, are in congruence with the computed data in the pDawn system. This leads us to either suspect, that the dark conditions in the pDusk system did not function or to go back to the model. We could conduct a parameter sensitivity analysis to find the most likely cause of the mismatch to data. With these results, we can nonetheless start and implement the equations for the mazEF mediated kill switch as we achieved a good fit to the more robust data given by Ohlendorf et al.¹.

Light on, BeeT off - Modeling The Optogenetic Kill Switch

After building the mathematical system for pDusk and pDawn, we were now able to simulate dark and light conditions. We have an understanding of how the optogenetic switch functions and can now explore the parameter space to find parameter sets which describe the wanted system behaviour including the toxin-antitoxin system. As there was no data available in literature to score the combined system, we used conservative constraints which are known from literature^{5 6} to construct our model for the optogenetic kill switch.

Goals

• Parameter Estimation: Find sets of parameters which give a response within the conservative constraints and thus show us where in parameter space to search in further analysis (e.g. relative kinetic rates of promoters).
• System Selection: %!here As the Biobrick ??? is not only used in pDawn, but also in the quorum sensing system, we wanted to have a tool to decide whether we can also use the less favourableBased on data from Ohlendorf et al. (2012), pDusk shows a smaller change in RFP_p ratio than pDawn. For the kill switch this translates into a smaller change in antitoxin over toxin (A/T) ratio and thus less likely to kill the cell effectively., but theoretically possible pDusk system in the kill switch design.

System Design

To built the kill switch we first studied the mazEF system and its functionality. The mazEF module is an antitoxin-toxin (TA) system native in E. coli. For the purpose of building a mathematical model, we describe the relevant aspects of the system here. For more background on the mazEF system, please go to the 5. The functionality of the system in nature is given by a fast degrading antitoxin and a slower degrading toxin. Once expression of the TA module is terminated, mazE will degrade faster and mazF will start cleaving mRNAs which ultimately leads to cell death. For the kill switch, we will consider either “under”-production of mazE in pDusk or overproduction of mazF in pDawn both under light inducing conditions (see Figure XY).

Constraints

For the parameter estimation of the optogenetic kill switch, all the assumptions mentioned in the above sections hold. Further, we assumed a lumped mazEF complex formation, where we disregard potential intermediary stages of complex formation. In addition, we defined the following new constraints:

• as we have literature data supporting that mazF degrades faster than mazE⁵, we exploited this as a constraint in our parameter generation procedure.
• another constraint was added to the system as the ratio of free A/T is assumed to be bigger than one⁶ under non inducing conditions.
• lastly, we added a logical constraint. Under the highest simulated light intensity, the A/T ratio has to be smaller than one

Results

We can show, that our parameter estimation procedure found two parameter sets, out of 1000 sampled sets, which satisfy the conservative constraints (see animated Figure 7 and Figure 8).

The two parameter sets give us an indication where in parameter space we will have to focus on for further analysis. These two sets also both describe a functioning pDusk and pDawn based optogenetic kill switch. Although, one might argue, that set number one is to be favoured as the change in mazE in pDusk compared to the change of mazE of set two in pDusk is higher and thus more likely to function and kill the cell. Future experiments on promoter strengths can be normalized to the model data. With such additional information, an indication can be given on which relative kinetic parameters are needed to make the system work.

When plotting the A/T ratio against the increasing light intensities, it becomes more evident, that set number one has a higher change in the A/T ratio within pDusk and is thus to be slightly favoured. To further evaluate these findings, more data from the lab is needed.

Methods

Looking at Figure_9 we can see, that adding the mazEF module to the pDusk/ pDawn system adds another 10 parameters. We constructed our parameter estimation similar to the parameter estimation of the previously described procedure and evaluated the system based on the constraints described above. It should be noticed here, that this is a qualitative evaluation and is not backed by statistical real world data. Nonetheless, it gives and indication about the real world situation.
We then took the best parameter set from the pDusk/ pDawn evaluation and added the new parameters to the ODEs. As we noticed that we have to discard around half of the parameter sets as they do not satisfy the constraint that mazE degrades faster than mazF, we built this constraint already into the parameter sampling procedure. In addition, a mazEF deficient E. coli. strain was engineered in the lab to prevent any interaction between the native and the engineered mazEF system.

%!here - missing ODEs Table of Parameter Set + Unit + explanation

Conclusion

%!here

System selection (cI_lambda inverter): As the Biobrick ??? is not only used in pDawn, but also in the quorum sensing system, we could now with additional lab data on different inverter strengths suggest which of those inverters will give the wanted response to balance the mazE mazF system.

Metabolic Modeling

In order to assess the real world viability of BeeT, we evaluated the proposed system of application by making a model of the entire system. To do this we used Flux Balance Analysis (FBA) to model the chassis. The chassis The chassis-organism is the framework that is modified for use in synthetic biology experiments. of BeeT is a variant of ​ Escherichia coli, for which it is known that it does not grow overnight in high osmotic pressure environments of 1 mol sucrose / liter or above. ¹ Supplementing an Apis mellifera (honey bee) colony with sugar water is a well established practice amongst beekeepers. ⁴ This is usually done with about 1 kg of table sugar for each liter of table sugarChemically speaking this is pure sucrose: C₁₂H₂₂O₁₁. After heating and stirring, this ends up as a concentration of 625 grams of sucrose per liter of water. This concentration can also be defined as 1.82 mol sucrose per liter, which is almost twice the threshold value at which E. coli would no longer grow. However, in this project we only need the ​E. coli to survive in the sugar water. This only needs to be as long as the worker bees require to pick the BeeT-cells up and bring them to where the Varroa destructor are. The question we try to answer is: does our chassis survive in the sugar water, and if so, for how long?

A model of E. coli survival

In order to make this model a technique called Flux Balance Analysis (FBA) was used to describe the relationship between the metabolism of the E. coli and the osmotic pressure of the sugar water. FBA is a mathematical method that can be used for predicting reaction fluxes and optimal growth rates of species using genome-scale reconstructions of their metabolic networks. In this case we used it to specifically predict the effect of water efflux on the total amount of ATP the cell has available for maintenance. The model used as a base for this analysis was iJO1366 ³from the BIGG database.

The relationship between maximum ATP available for survival and water efflux is shown in Figure 1 the results from the metabolic model suggest that there is a linear relationship between ATP availability and water efflux. This implies that if no water is available for ATP used for maintenance, the cell will die. When the model runs without any modification, i.e. in an environment where it is in the exponential growth phase, an ATP Maintenance flux of 3.15 mmol⋅gram Dry Weight^-1⋅hour^-1 is given as output by the model.

We do not know the exact amount of ATP needed for survival in sugar water conditions, but because of the results from Figure 1 we can start looking at the relationship between survival time and water efflux.

In Figure 2 we can see not only the relationship of survival time against max ATP available for survival, but also how different thresholds of minimal cell-water tolerance would affect this relationship. The minimal cell-water tolerance threshold gives the value at which percentage of the remaining cell-water reaches a state that will cause cell-death. This has a drastic effect on survival time, changing 20 minutes of maximum survival time to a mere ~90 seconds in the worst case scenario.

Up to 90 minutes

Figure 1 shows us that osmotic pressure alone can indeed have an effect on cell regulation and cell death, and from Figure 2 it appears that the minimal water allowance threshold has a high impact on range of possible times. Our model is limited in that it predicts an infinite survival time beyond 90 minutes. This suggests that our model may be missing some form of regulation that allows for longer survival times. To understand this effect, we conducted an experiment to test this prediction. This experiment can be found here.

What we can say is that if the cells can survive for longer outside of this period, then they must have enough ATP available for basic maintenance. If cell death occurs, other processes than pure water-efflux must be the cause of that. This could be due to a combination of lacking nutrients and water-efflux, or over production of osmolytes to keep the balance.

Beehave

Due to regulatory and experimental hurdles it is difficult to test the effectiveness of BeeT in combating Varroa destructor in the field. We would still like to be able to give advice to local beekeepers on the ideal application strategy of BeeT based on several scenarios. To accomplish this, the well-known model BEEHAVE was adapted to include the effect of BeeT on Varroa mite and virus dynamics in simulated colonies. BEEHAVE is an open-source, agent-based model which can be used to examine the multifactorial causes of Colony Collapse Disorder ^{1 2} . It consists of several modules which controls aspects like foraging, Varroa mite dynamics and colony growth (figure 1) and was extensively tested for robustness and realism ¹ .

BeeT can be transported into the hive by adding it to sugar water or bee bread. Each of these sources has its own advantages and disadvantages. In this project we have used Escherichia coli as our model chassis and this can be applied to sugar water [link to Ronald]. Sugar water is supplemented to support a colony during honey harvest and before winter as a substitute for nectar ^{3 4} . Supplementing Apis mellifera (honey bee) with sugar water is a well established and familiar practice amongst beekeepers ³ . However, reconstructing the same system in Lactobacillus species would allow us to use BeeT in bee bread. Bee bread has certain advantages over sugar water as it is more frequently transported to larvae. It is a combination of pollen, regurgitated nectar and glandular secretions and is inoculated with fermenting bacteria by honeybees ⁵ . There is mounting evidence that pollen supplementation increases protein content in honey bee hemolymph, likely improving survival of colonies to various stressors ^{6 7} .

BeeT Module

BEEHAVE is an open-source, GNU licensed agent-based model utilizing NetLogo and consists of several interlocking modules which each model different aspects of the beehive. BEEHAVE has the intended goal of modelling the wide variety of stresses affecting honey bees and is the only model incorporating all these different aspects ¹ . As such it is the ideal basis for our investigation into the effects of BeeT on Varroa mites and honey bee dynamics. BEEHAVE has several modules covering inter-colony dynamics, foraging and a Varroa mite model as depicted in Figure 1. Two viruses are also included in the model; deformed wing virus (DMV) and acute paralysis virus (APV) for which Varroa mites are a vector. Our BeeT module, which runs parallel to BEEHAVE, is capable of modeling transport of BeeT into the hive using sugar water or bee bread. It also calculates how much BeeT is transported to larvae based on consumption of respectively honey and pollen stores. Based on the amount of BeeT at larvae the Varroa mite mortality, when Varroa mites emerge from brood cells, is determined. This in turn affects Varroa mite population levels in the hive, reducing virus loads in the hive and allowing colony survival.

Assumptions

Not everything is known about the BeeT model and as such it was necessary to make certain assumptions. Some of these assumptions are related to BeeT in general, while others are specific to the sugar water or bee bread applications. We will first discuss general limitations of our model and then address each of the applications separately.
General assumptions
Several properties of BeeT and its chassis are currently unknown. This includes the degradation rate of BeeT both inside and outside the hive. With both treatments BeeT encounters different environments while being transported to larvae. For honey, BeeT is exposed to sugar water, the sugar stomach of bees, honey stores and brood cells. While bee bread encounters pollen and brood cells in its travels. We do not know the effects of these different environments on degradation rates of BeeT. For this reason we assumed that degradation rates are stable inside and outside the hive.
Another uncertainty is what the effect of BeeT is on Varroa mite mortality, we modelled this by using a saturating function. We also assumed that the effect of BeeT on Varroa mite mortality is entirely determined by the amount of BeeT at the larvae at the time of capping a brood cell.
Finally, the BEEHAVE model is only able to model a single virus, either DMV or APV. Consequently, we are unable to model the combined effects of both viruses on a bee colony. For all analysis, we used the DWV virus, as it is more harmful to honey bee colony survival than APV 8 .
Sugar water
We know that, based on the experimental results [hyperlink] obtained by Ronald, E. coli is able to survive in sugar water for at least 48 hours. His model additionally predicts that if E. coli is capable of surviving 48 hours, it can likely survive long enough to be transported into the hive [hyperlink]. His results do not predict the degradation rate of E. coli in sugar water, honey or in brood cells.
Additionally, we estimate that in sugar water there are approximately 1 * 10^6 cells*ml^-1. This estimate is based on the assumption that we add saturated medium (1 * 10 ^ 8 cells*ml^-1) containing BeeT to sugar water and diluted this by a factor of 100. The dilution is performed to avoid rejection of sugar water by bees as it would be too contaminated to consume.
We based the hourly uptake rate of sugar water on a paper by Avni et al. ⁹ , it is possible that this is not the maximal uptake rate of sugar water. In the measured period all of the sugar water was taken up by the colony, making a higher uptake rate possible.
Bee bread
We assume that the amount of cells in bee bread is roughly equal to the number of CFU*ml^-1 in yoghurt ¹⁰ . Artificial bee breads can be made using yoghurt and these are readily accepted by bees ¹⁰ . Additionally, we assume that the removal rate of artificial bee bread is equal to 22.7 g*day^{-1 4} . Based on Avni et al. we know that the uptake rate depends on the manner in which bee bread is applied ⁹ . We chose to base the uptake rate of bee bread on Brodschneider et al. since their experimental setup is similar to Avni et al. ^{4 10}

Key results

The functionality of the BeeT module is primarily based on literature research and as such requires a relatively low amount of parametrization. The main unknown quality is how BeeT will behave, namely degradation and its effectiveness at combating Varroa mites at larvae. We are interested in how BeeT can best be applied given certain assumptions [hyperlink previous section]. As such we are interested in the best period to apply BeeT; either before winter or during spring so we can give recommendations to beekeepers. We also examine the difference in effectiveness between BeeT in sugar water and in bee bread for both periods. This can inform future work on whether it is beneficial to adapt BeeT for functionality in bee bread. There are three additional parameters, related to BeeT, that can be varied and upon which effectiveness relies. These are; degradation of BeeT in the hive, outside the hive and the effect of BeeT on Varroa mite mortality and each is varied across a range of values. Additionally, for every combination we have 5 replicates to reduce the variance in our results. We divided the analyses into 4 treatments, each with the same three parameter ranges and 5 replicates per set of parameters. This can be seen in table 1.

Period	Treatment	Range
1st April - 15th June	Sugar Water	Varroa Mite Mortality, Degradation in-hive and Degradation outside-hive
	Bee Bread	Varroa Mite Mortality, Degradation in-hive and Degradation outside-hive
1st September - 15 November	Sugar Water	Varroa Mite Mortality, Degradation in-hive and Degradation outside-hive
	Bee Bread	Varroa Mite Mortality, Degradation in-hive and Degradation outside-hive

Each combination of period and treatment has the same ranges. We compared the winter period (1 September - 15 November) with the spring period (1ste April - 15 June) for the different ranges and treatments. In each case the spring period was significantly more effective in combating Varroa mites than the winter period. For further analysis we continued with the spring period.

Three scenarios were identified; the colony dies due to virus load, the colony survives at a lower population level or the colony thrives. Every year the BEEHAVE model checks the population after overwintering, if it is below 4000 worker bees the colony is presumed dead. To differentiate between the different scenarios we looked at the mean population of worker bees after overwintering. If this number is below 3000 mean overwintering bees we assume that the colony is dead. If the mean population of overwintering bees is between 3000-5000, the colony is capable of surviving. Above 5000 overwintering worker bees, the colony is thriving and the Varroa mite infestation is under control.Each parameter combination, for sugar water and bee bread, falls within one of these three scenarios. For both treatments, we chose a parameter set which is representative for a certain scenario. From this we have 6 representative parameter sets; three for sugar water (death,survival and thriving) and three for bee bread (death,survival and thriving). As can be seen in table 2, bee bread is significantly more effective than sugar water. Some BeeT is consumed by worker bees instead of deposited at the larvae. This is represented by the fraction of total BeeT moved to the hive that is consumed by larvae and also depicted in table 2.

In the case of sugar water a far lower fraction of BeeT ends up at larvae and most is consumed by worker bees. This is likely the main reason for the difference in effectiveness for sugar water and bee bread.

Each of the 6 parameter sets was further analysed by examining the population dynamics of worker bees and Varroa mites over time. In figure 2 and figure 3 these are shown, respectively, for a colony barely surviving a Varro mite infestation and thriving regardless of the presence of Varroa mites at the start of the simulations. When the colony is barely surviving we can see that worker bee and Varroa mite populations reach a dynamic equilibrium.

When the colony is thriving, the starting Varroa mite population rapidly declines and disappears.

This indicates that the treatment is effective for a low population of initial Varroa mites, we further analysed a situation with a higher starting population of Varroa mites which is shown in figure 4.

Even with a higher starting population and a high initial virus load, we can see that the treatment is successful in reducing Varroa mite populations and that the colony recovers.

Model analysis

The BEEHAVE model utilizes multiple pseudo-random number generators to calculate mortality amongst Varroa mites and bees. These pseudo-random number generators increase the variance inherent in the model and as such it is advisable to use multiple replicates per parameter set. To determine the error bars we ran a simulation where we varied the Varroa mite mortality and plot the mean worker bee population over 10 years with 10 replicates per value of Varroa mite mortality. This was done for the first version of the BEEHAVE model, which governs year round Varroa mite mortality when they emerge from a brood cell.

With 10 replicates the variance is relatively low, particularly when the mean worker bee population is high.
To analyse the full BeeT module we chose to use 5 replicates, this is a compromise between simulation time (137.100 simulations with 5 replicates running for ~ 8 days) and reduction of variance.

We first separated the results of the analyses of the full BeeT module based on the three scenarios. This resulted in figure 6, this is a 3d graph with the axis degradation in the hive, degradation outside the hive and impact of BeeT on Varroa mite mortality. Additionally, the size of each set of parameters indicates the mean overwintering population of worker bees.

Discussion and Conclusions

Currently it is unclear how effective BeeT will be at combating Varroa mites, primarily toxicity is difficult to measure in the lab and deploying BeeT in the field has many regulatory issues. As such we build a BeeT module for the well known BEEHAVE model to assess the minimum effectiveness for certain treatments. We also examined two different periods for administering the BeeT; before winter and during spring. Before winter was chosen to reduce the Varroa mite population before winter and support honey bees during their most fragile period in the year (winter). Additionally, it is normal for current beekeeper practices to supplement with sugar water during this period and this period would provide a minimal amount of disturbance for beekeepers. The advantage of administering BeeT during spring is that this is the period at which most of the Varroa mite population is present at larvae. By administering during this period we can ensure that BeeT is able to affect a large percentage of the Varroa mite population. Beekeepers also visit hives regularly during this period as they harvest honey in spring. Recent papers indicate that feeding of artificial bee bread during this period can have significant positive effects on honey bee health ^{6 7} . For both sugar water and bee bread the ideal period for administering BeeT is during the spring, due to the large difference in effectiveness between the two periods. In a future application of BeeT this can be given as a recommendation to beekeepers. This is a practice which can be easily adopted since beekeepers already visit hives regularly in the spring period. In particular for supplementation of sugar water, as this is already accepted practice amongst beekeepers.

The main difference between sugar water and bee bread is that a larger fraction of total in-hive BeeT is given to larvae for bee bread. This is not surprising as honey stores are primarily used to support energetic needs of worker bees ¹ and as such only a small fraction is transported to larvae. On the other hand, the protein contained in bee bread is crucial for the proper development of larvae and as such is mainly given to larvae to support their growth ¹ . Based on our results we expect that BeeT will be effective using sugar water but greater effectiveness can be accomplished by engineering BeeT for compatibility with bee bread.

Modeling methods

In version one of the model the following code was added to BEEHAVE:

This calculates the mortality of Varroa mites when they emerge from the brood cell (mortalityMites). In BEEHAVE the mortality of Varroa mites is encapsulated in the fecundity of female Varroa mites, this means that mortalityNoBeeT is zero. When ScalingVariable is large, the mortalityMites approaches one.
The amount of Varroa mites emerging from a cell is calculated with the variable ‘healthyMitesInSingleCell’, using a random-poisson distribution some of the Varroa mites emerging from the brood cell die. Both versions of the BeeT module only affect the Varroa mite population when Varroa mites emerge from a brood cell. This was done to reduce the chance of introducing bugs into the original BEEHAVE model.

In version 2 of the BeeT module, the mortalityMites was changed to:

So that it is dependent on the concentration of BeeT in the cell at the time of capping. Additionally we included a procedure “BeeTProc which is called daily. It handles BeeT moving from outside to inside the hive at a set rate per hour. It also handles degradation, both hourly degradation outside the hive and daily degradation inside the hive. When Varroa mites invade, the amount of BeeT present at larvae is saved for each invading Varroa mite in the variable BeeTLarvaeInvasion. Whenever a Varroa mite emerges from a cell its mortality is calculated based on BeeTcell which is dependant on BeeTLarvaeInvasion. The total amount of BeeT moving to larvae is based on the amount of sugar water and pollen consumed by larvae. If BeeT is consumed by worker bees it is considered lost and is removed from the BeeT in-hive store.

All variables added to BEEHAVE can be found in table 3, including how they are calculated, if applicable what their initial value is, units, in what procedure they are used and any comments or references.

References

1. Eri Nasuno, Nobutada Kimura, Masaki J. Fujita, Cindy H. Nakatsu, Yoichi Kamagata, and Satoshi Hanada (2012). Phylogenetically Novel LuxI/LuxR-Type Quorum Sensing Systems Isolated Using a Metagenomic Approach Vol, 78, number 22. ↩

2. Author unknown, Retrieved from ↩

3. Pearson, B. Bacterial Hash Function Using DNA-Based XOR Logic Reveals Unexpected Behavior of the LuxR Promoter, IBC 2011, vol. 3, article no. 10, pp. 1-10 ↩

4.Sergi Abadal Ian F. Akyildiz (2011). Automata modeling of Quorum Sensing for nanocommunication networks, Volume 2, Issue 1, Pages 74–83 ↩

5.Designed by: Kin Lau Group: iGEM09_MoWestern_Davidson (2009-07-22) ↩

6. Gregory Bokinsky, Edward E. K. Baidoo, Swetha Akella, Helcio Burd, Daniel Weaver, Jorge Alonso-Gutierrez, Héctor García-Martín, Taek Soon Lee, and Jay D. Keasling(2011).HipA-Triggered Growth Arrest and β-Lactam Tolerance in Escherichia coli Are Mediated by RelA-Dependent ppGpp Synthesis 2013 Jul; 195(14): 3173–3182. ↩

7. A. Raue et al. Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Vol. 25 no. 15 2009, pages 1923–1929 ↩

8. Martin Bengtsson et al. Gene expression profiling in single cells from the pancreatic islets of Langerhans reveals lognormal distribution of mRNA levels volume 15 page 1388–1392(2005) ↩

1. Cheng, Y. L., Hwang, J., & Liu, L. (2011). The Effect of Sucrose-induced Osmotic Stress on the Intracellular Level of cAMP in Escherichia coli using Lac Operon as an Indicator. Journal of Experimental Microbiology and Immunology (JEMI) Vol, 15, 15-21. ↩

2. Neidhardt F.C. Escherichia coli and Salmonella: Cellular and Molecular Biology. Vol 1. pp. 15, ASM Press 1996 ↩

3. Orth, J. D., Conrad, T. M., Na, J., Lerman, J. A., Nam, H., Feist, A. M., & Palsson, B. Ø. (2011). A comprehensive genome‐scale reconstruction of Escherichia coli metabolism—2011. Molecular systems biology, 7(1), 535. ↩

4. DeGrandi-Hoffman, G., & Hagler, J. (2000). The flow of incoming nectar through a honey bee (Apis mellifera L.) colony as revealed by a protein marker. Insectes Sociaux, 47(4), 302-306. ↩

1. M. A. Becher, V. Grimm, P. Thorbek, J. Horn, P. J. Kennedy, and J. L. Osborne, (2014). BEEHAVE: A systems model of honeybee colony dynamics and foraging to explore multifactorial causes of colony failure, Journal of Applied Ecology., vol. 51, no. 2, 470–482. ↩

2. J. Horn, M. A. Becher, P. J. Kennedy, J. L. Osborne, and V. Grimm (2015), “Multiple stressors: Using the honeybee model BEEHAVE to explore how spatial and temporal forage stress affects colony resilience,” Oikos, no. September 2015, 1001–1016. ↩

3. Fasasi, K A (2011) "Cumulative effect of sugar syrup on colony size of honeybees, Apis mellifera adansonii (Hymenoptera : apidae) in artificial beehives " Journal of Natural Sciences, Engineering and Technology ed. 10: 33-43. ↩

4. Robert Brodschneider, Karl Crailsheim (2010) "Nutrition and health in honey bees" Apidologie ed. 41: 278-294 ↩

5. Ellis, Amanda M. Hayes, Jr, G. W. (2009) "An evaluation of fresh versus fermented diets for honey bees (Apis mellifera)." Journal of Apicultural Research 48: 215-216 ↩

6. Basualdo, Marina Barragán, Sergio Antúnez, Karina (2014) "Bee bread increases honeybee haemolymph protein and promote better survival despite of causing higher Nosema ceranae abundance in honeybees" Environmental Microbiology Reports 6: 396-400 ↩

7. De Jong, David (2009) "Pollen substitutes increase honey bee haemolymph protein levels as much as or more than does pollen" Journal of Apicultural Research Reports 48: 34-37 ↩

8. Sumpter, D. J T, Martin, S. J., (2004). The dynamics of virus epidemics in Varroa-infested honey bee colonies, Journal of Animal Ecology., vol. 73, is. 1, 51-63. ↩

9. Sumpter, D. J T, Martin, S. J., (2004). The effect of surface area of pollen patties fed to honey bee (Apis mellifera) colonies on their consumption, brood production and honey yields., Journal of Apicultural Research., vol. 48, is. 1, 23-28. ↩

10. Ellis, Amanda M. Hayes, Jr, G. W., (2009). An evaluation of fresh versus fermented diets for honey bees (Apis mellifera)., Journal of Apicultural Research., vol. 48, is. 3, 215-216. ↩




iGEM Wageningen UR 2016. Send us an email or follow us on Facebook!
For more information about iGEM, check out their website and Twitter!