Team:Wageningen UR/Model

MENU ▤
Wageningen UR iGEM 2016

 

Modeling Overview

In light of our guiding principles (specificity, regulation and biocontainment), we modeled four different aspects of BeeT. The modeling work can inform and aid the improvement of wet-lab experiments. Another facet is to assess the optimal application strategy for BeeT in real world applications. We asked ourselves; what parts of our system can benefit the most from an interplay between modeling and experimental work? These considerations led us to ask the following questions:

Population Dynamics

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 modeling.

Quorum Sensing Circuit

Quorum sensing is a cell-cell communication system. The detection of chemical molecules allows the bacteria to sense the presence of other bacteria. In this way the bacteria control gene expression in response to changes in cell number1,4. This process is achieved through the production and release of an Acylated Homoserine Lactone (AHL) autoinducer An autoinducer is a signal molecule related to cell density. . The AHL can diffuse 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.
The quorum sensing system consists of two proteins: LuxI and LuxR. LuxR is constitutively expressed, together with AHL it forms a complex which activates the transcription of a toxin. In our case we use the detection protein Green Fluorescence Protein (GFP) to follow the behavior of the system2. LuxI encodes an AHL synthase. AHL is a molecule that can diffuse freely through the cell membrane, and in this way travels from cell to cell. LuxR encodes for a protein and together they form a complex. In a natural system the LuxR-AHL complex controls transcription of LuxI, forming a positive feedback loop that increases the amount of AHL in the system3. When there is more AHL in the system, AHL is more likely to bind to the LuxR protein. In Figure 1 you can see a schematic representation of the system.

Figure 1: Schematic of the quorum sensing system. The arrows indicate activation of a part and the flat bars indicate inhibition of a part. The quorum sensing system consists of two parts: LuxR and LuxI. LuxR constitutively produces a transcriptional factor that binds to AHL and LuxI produces an AHL autoinducer that diffuses freely through the cell membrane, from cell to cell.

According to the 2011 team from Davidson College and Missouri Western State University5, LuxR protein represses transcription of the luxR gene in the absence of AHL, establishing a negative feedback. We investigated whether this negative feedback can be used to create different populations. As you can see in Figure 3 different states of GFP responses can occur when the production rate of luxR and rate of complex formation are changed. However, when we remove the feedback in the system and change the same parameters, we get similar responses of the system as shown in Figure 2.

Figure 2: Quorum sensing without negative feedback. In the left graphs you see the production of GFP when the production rate of LuxR and complex formation parameters are different in each cell. Every line represents one cell, the system is simulated with 10 cells. Five of the 10 cells have a value of 0.01 for the production rate of LuxR and the complex formation parameters, and the other half of the cells have a value of 1.0 for these parameters. In the right graphs you see the responses of the GFP where the parameter values are the same for each cell.
Figure 3: Quorum sensing with negative feedback. In the left graphs you see the production of GFP when the production rate of LuxR and complex formation parameters are different in each cell. Every line represents one cell, the system is simulated with 10 cells. Five of the 10 cells have a value of 0.01 for the production rate of LuxR and the complex formation parameters, and the other half of the cells have a value of 1.0 for these parameters. In the right graphs you see the responses of the GFP where the parameter values are the same for each cell.

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

Figure 4. In this graph you see the amounts of GFP we got from the different parameter sets for the system without negative feedback. Most of the parameter sets give a low GFP response, but a few give high ones. Those we use for further analysis.

Figure 5. In this graph you see the amounts of GFP we got from the different parameter sets for the system with negative feedback. Most of the parameter sets give a low GFP response, but a few give high ones. Those we use for further analysis.

The figures are simulated with the parameter sets obtained from the quorum sensing model. We are using the parameter sets that give the best response according to the confidence intervals (see Method section). In the presence and absence of negative feedback, the system behaves similarly by assuming the production rate of luxR and the formation rate of the complex are different in each cell. The quorum sensing system can simulate different responses with the same cells and different parameters, if we could control the LuxR production rate. However, we could not envision a method of doing this experimentally. We aim to make cells that are genetically identical but still can respond differently. Therefore, we designed another system that could help us create this response downstream of the quorum sensing mechanism. The two systems were tested in the lab parallel with the modeling.

Why Include a Subpopulation System?

Since the quorum sensing mechanism could not provide the subpopulations we desired using genetically identical cells, we created a system that could act downstream of quorum signalling to provide the response we needed. Population-wide Cry toxin overexpression is likely to kill all E. coli cells in the first moment of toxin expression. Since this would only affect one of the subpopulations, the second subpopulation would be able to initiate a new growth phase after death of the toxin-producing cells. The critical requirement for this is that cells respond at different times to the quorum stimuli despite being genetically identical. A small part of the population that acts differently from the rest of the population is called a subpopulation. 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 cI inactivates the λ cI directly, or prevents the translation of the λ cI protein. 434 cI has a higher turnover rate than λ cI. This subpopulation system is based on a natural system for persister cell formation 6.

Figure 6: Schematic of the subpopulation system. In this system, AraC inhibits the promoter activity of the subpopulation part. AraC is inhibited by Arabinose, and when glucose is added to the system Arabinose is inhibited by glucose. AraC inhibits the constitutive promoter of 434 cI and λ cI, where 434 cI inhibits the promotion of λ cI. The production of λ cI leads to production of Red Fluorescence Protein(RFP) as the reporter protein.

As you can see in Figure 6; glucose has a suppressing function on the system, arabinose has an 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 cI degrades faster than λ cI, cell division happens when the cell size has doubled, and two cell populations exist with one growing differently than 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. To see how the model behaved when we changed the ratios of the initial amounts of λ cI and 434 cI, we investigated the effect of starting conditions on the outcome of the system. The initial conditions were varied between 1 and 10. Within the heatmap you can see in which ratios the initial amounts of λ cI and 434 cI in the system are needed to get high RFP production (here glucose and arabinose have a fixed concentration).

Figure 7. Heatmap of λ cI and 434 cI against RFP production. The right side bar indicates the amount of RFP, where red indicate a high amount and black a low amount. This heatmap shows you the amount of RFP produced given the balance between the two genes initial conditions.

As seen in Figure 7, when λ cI and 434 cI are present in similar amounts we have a subpopulation system. But if they are present in different amounts there will be no subpopulations. This can be expected when you look at the subpopulation system. The system is inhibited by 434 cI, which represses the RFP production, and λ cI activates the RFP production. In Figure 7 you can see the balance between the 434 cI and λ cI amounts that are present for the output of RFP. To understand how changing the Ribosomal Binding Site (RBS) impacts RFP responses, we simulated the system for different combinations of transcription rates. From these simulations we found:

Table 1: RBS library of the 434 cI gene with different predicted translation rates.

λ cI 434 cI RFP responses
1.94 135.6 7.00
1.94 108.8 8.51
1.94 942.5 1.48
1.94 726 1.77
1.94 628.6 1.962
1.94 577.1 2.09
1.94 529 2.23
1.94 40.3 20.1
1.94 384.9 2.87
1.94 268.5 3.87
1.94 257.8 4.01
1.94 256.7 4.02
1.94 73.2 12.1
1.94 100.9 9.09

To confirm that our model represents biological behavior we analysed the RBS library. We can conclude from these results that the RFP concentration qualitatively corresponds to the Figure 7 Because as you can see in the Table 1 the highest RFP responses are when the λ cI and 434 cI are present in the same amount. The lowest responses of RFP when there is a is a big difference between λ cI and 434 cI. This is also shown in the heatmap we got form the model.

Combined System

Figure 8: Schematic of the quorum sensing system attached to the subpopulation network. The arrows indicate activation of a part and the flat bars indicate inhibition of a part. The quorum sensing system consists of two parts: LuxR constitutively produces a transcriptional factor that binds to AHL, and LuxI produces an AHL autoinducer that diffuses freely through the cell membrane, from cell to cell. LuxR forms a complex with AHL that inhibits the subpopulation system, where 434-cl-LVA and λ cl are under the same promoter. This leads to either GFP production, or RFP production depending on the same strength of quorum sensing.

Based on the modeling we hypothesized the following: when there are more cells present in the system, more AHL-LuxR complex is formed. The complex inhibits the subpopulation promoter. When the promoter is inhibited production of 434 cI will be suppressed and production of λ cl will be activated. At a certain time point λ cl takes control over the system, because 434 cI 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 intend to generate different cell populations. With this extended network we try to generate a model which predicts the behavior of the subpopulations. In Figure 9 you can see the increasing RFP over time after many cell divisions, indicating an increasing cell population.

Figure 9. Volume oscillations correspond to the dividing cells. Together they form a growing cell population. On the left y-axis you see the RFP, and on the right y-axis you see the volume. The x-axis is the time. The total amount of RFP produced by the toxin is shown by the red line.

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 between zero and one. For the analyses we used 100,000 parameter sets, 10 cells, random initial conditions between zero and one 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 HypercubeLatin hypercube sampling is a statistical method to sample the parameter space equally. It divides the parameter space according to the number of parameters and sample size and picks a random number from each parameter subspace. sampling.

With 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 distribution8 with a parameter estimation based on Raue et al7. In Figure 10 you can see an example.

Figure 10. Here you see the lognormal distribution of the 19 parameters in the quorum sensing system. The red line indicates the mean of the parameter.

With these lognormal distributions and the accessory confidence intervals the best parameter sets could be chosen.

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 behavior of the subpopulations. In figure 9 you can see the increasing RFP over time after many cell divisions, indicating an increasing cell population.

Click the button to see the equations!

Conclusion

Quorum sensing ensures that the toxin is only produced when the density of bacteria is high enough, thus standardising the amount of toxin produced by the bacteria population. The subpopulation system was coupled to the quorum sensing system. Together, quorum sensing and formation of subpopulations allow bacteria to produce waves of toxin.
Different parameters in the quorum sensing system elicit different GFP responses. However, we can not expect to be able to make genetically identical cells with varying values for these parameters With the subpopulation system we may be able to predict which RBSs of the library have the highest chance to successfully create subpopulations. The library can be tested and used in further research.


Optogenetic Kill Switch Model

One important aspect of our project is the biocontainment of BeeT to the hive, preventing it from spreading into the environment. In collaboration with the wet-lab, we designed an optogenetic kill switch based on the blue light sensing pDusk and pDawn systems as described by Ohlendorf et al1. The cell death is inflicted by the mazEF toxin-antitoxin (TA) system. The latter functions by the antitoxin (mazE) preventing the toxin (mazF) from cleaving mRNAs which ultimately kills the cell. With the underlying biology of the system in mind, we 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 al1. 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 tools.

Light On

Modeling The Optogenetic Tool

For the first part of our analysis, we optimized the mathematical models of pDusk and pDawn using published data1. 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.
By constructing an accurate model of the optogenetic part of our overall kill switch system, we are confident that we are mathematically representing the behavior of the microorganism in the lab. Results from the model can then be used as a reference point for building the toxin-antitoxin system into the kill switch model. We will search the parameter space to identify parameter sets which fit the data given by Ohlendorf et al1. In addition, we will select the best fit to data based on a scoring function as described by Raue et al4. Lastly, we provide a mathematical model which can be used by future iGEM teams to build new optogenetic tools.

System Design

We first studied the pDusk system in more detail to derive simplified equations which can describe the system behavior. The genetic switches are based on the interactions between the light sensing protein YF1 and the response regulator FixJ as described by Möglich et al2. In darkness YF1 phosphorylates FixJ, which in turn activates the promoter pFixK2. This activation leads to the expression of the 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 pFixK2 promoter and thus represses the target genes in Figure 11 for pDusk.

Figure 11. Simplified representation of the pDusk and pDawn system with an RFP reporter gene. The arrows indicate activation of a part and the flat bars indicate inhibition of a part.

The pDawn system functions in the same way as the pDusk system, with the addition of λ cI being expressed depending on the activity of the pFixK2 promoter. The protein form of λ cI inhibits the expression of the RFP reporter gene in Figure 11 for pDawn.

Assumptions

To derive mathematical equations describing the above systems, we made some simplifying assumptions.

  • We neglect the effects of dilution/ growth. We thus assume the cell’s volume to be constant over time.

  • We assumed that the components in our system are neither affected by nor affect other cellular mechanisms.

  • As there is no additional information gained from modeling transcription and translation for the dark state of YF1(yDD) and the inactive form of FixJ (ji), we lumped these processes and assumed that they can be described by a single parameter.

  • We also assumed the rates of production, binding, and degradation of mRNAs and proteins to be linear which results in a first order differential equations system

To model the influence of light on the system, we assumed light activation to occur similar to the model of phytochrome B dimerization described in Klose et al3 as N⋅k 2. In addition, we consider in our system the different stages of YF1 (yDD, yDL/LD, and yLL) as described by Möglich et al2. Translated to our system, this means that with increasing light intensities (N), the pool of activated YF1 (yDL/LD and yLL) would deplete the dark state of YF1 (yDD), which is responsible for the phosphorylation and thus activation from inactive FixJ (ji) to active FixJ (ja). As all stages of YF1 form a dimer, we assumed dimerization to occur quickly 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 (ja) on the expression downstream of the RFP (RFPm) reporter gene. A detailed system design is provided here.

Results

With our mathematical model we can describe the behavior, as shown in Figure 12, of both systems as tested in the lab by Ohlendorf et al1. We can thus say that we built a simplified mathematical representation of the pDusk and pDawn system which can describe the systems’ behavior.

Figure 12. Comparison of computed model results to the data of Ohlendorf et al1. 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 al1. 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 behavior of the system.

Preliminary results from the lab showed that the RFPp response after 17 hours incubation time is weak using the pDusk system. The experiment was conducted in darkness and under 1.41 ⋅ 10 4 [ µ m o l ⋅ m - 2 ⋅ h - 1 ] light intensity and the results are shown in Table 2.

Table 2: Preliminary Lab Results on pDusk/ pDawn System
System Dark Light
pDusk 45.43 ± 37.83 [ F L O D 600 ] 36.79 ± 18.35 [ F L O D 600 ]
pDawn 126.54 ± 25.19 [ F L O D 600 ] 5356.75 ± 78.55 [ F L O D 600 ]

This data from the lab suggests that there is not a high change in RFP for the pDusk system between dark and light conditions. When plotting the computed response of RFP for increasing light intensity in Figure 13 we see, however, that pDusk has a steep decline in RFP levels from complete darkness to very low light intensities. There could be two potential reasons for such a result. First, despite all efforts in creating dark conditions, the experiment could have been contaminated with low light intensities. This would explain why there is no significant change visible in the experimental data between different conditions. To confirm this, additional experiments would need to be conducted. Secondly, a more in-depth parameter sensitivity and robustness analysis could also reveal potential causes of the discrepancies between model simulations and experiments conducted under darkness. This could lead to an understanding of how to engineer the promoters in pDusk to amplify the wanted the response.


Figure 13: RFP Response in pDusk and pDawn per Light Intensity.

Methods

After studying the pDusk and pDawn system, we constructed a parameter estimation procedure by picking random parameter sets based on the latin hypercubeLatin hypercube sampling is a statistical method to sample the parameter space equally. It divides the parameter space according to the number of parameters and sample size and picks a random number from each parameter subspace. sampling (LHS) method in Matlab 2015b (used for the whole study). We simulated 10,000 parameter sets and assumed the initial conditions of only yDD and ja being in the on state, i.e. we assigned the value one. All other initial conditions were kept at zero. 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 lognormal distribution that is representative of measured biological rates7. We transformed the two parameters, k2 and k3, for which experimental data was available according to their known distribution. We then loaded the RFP response data ± standard deviation from Ohlendorf et al.1 for both pDusk and pDawn and transformed the light intensities from [µW⋅cm-2] to [µmol ⋅m-2⋅h-1]. For a more detailed explanation see the notebook.

We then simulated both systems according to our ordinary differential equations (ODEs) with the same parameter sets, and stored the response at 17 hours for later evaluation in the scoring function. We simulated the system for 17 hours, as this is the duration Ohlendorf et al.1 used in their study. For the next step, we normalized the computed RFPp response with the measured response at the highest light intensity used by Ohlendorf et al1. A more detailed description of how we constructed the scoring function and the ODEs for both systems can be found here:

Click the button to see the detailed methods!

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 tool functions and can now explore the parameter space to find parameter sets which describe the wanted behavior and regulation of the toxin-antitoxin system. As there was no experimental data available in literature to score the combined system, we used conservative constraints which are known from literature5,6 to construct our model for the MazE and MazF kinetics.
The objective is to 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 to describe the toxin-antitoxin system. As the Biobrick BBa_K081007 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 sensitiveBased on data from Ohlendorf et al. (2012), pDusk shows a smaller change in RFPp ratio than pDawn. For the kill switch this translates into a smaller change in antitoxin over toxin (A/T) ratio and thus a smaller likelihood to kill the cell effectively., but theoretically feasible pDusk system in the kill switch design.

System Design

To build the kill switch we first studied the mazEF system and its functionality. The mazEF module is an toxin-antitoxin (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 wet-lab part.
In the TA system, mazE represents the antitoxin which forms a dimer and binds two dimers of mazF to built an inactive complex under non-lethal conditions5. 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 as shown in Figure 14.

Figure 14. Simplified representation of the pDusk + constitutive mazF and pDawn + constitutive mazE system. The arrows indicate activation of a part and the flat bars indicate inhibition of a part.

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 intermediate stages of complex formation. In addition, we defined the following constraints:

  • As we have literature data supporting that mazF degrades faster than mazE5, 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 one6 under non-inducing conditions.

  • Lastly, we added a logical constraint which describes the underlying biology of the system. Under the highest simulated light intensity, the A/T ratio has to be smaller than one to show that toxin wins over antitoxin and kills the cell.

Results

We can show that our parameter estimation procedure found two parameter sets, out of 1,000 sampled sets, which satisfy these constraints as shown in Figure 15 and 16.

Figure 15. Response of the two parameter sets which satisfy the conservative constraints. The legend in pDusk is also valid for pDawn.

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 (blue and red lines) is to be favoured as the change in mazE in pDusk compared to the change of mazE of set two (yellow and purple lines) in pDusk is higher and thus more likely to function and kill the cell. With this information, an indication can be given on which relative kinetic parameters are needed to make the system work.

Figure 16. A/T ratios per light intensity of the two parameter sets which satisfy the constraints. Where blue = set 1 and red = set 2.

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

Methods

Adding the mazEF module to the pDusk/ pDawn system requires another 10 parameters to describe the system. 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 noted, that this is a qualitative evaluation and is not backed by statistical real world data. Nonetheless, it can be used as an indication.
We then took the best parameter set from the pDusk/ pDawn evaluation and added the equations describing the interactions with the MazEF system to the set of ODEs. As we noticed that we have to discard around half of the parameter sets because they do not satisfy the constraint that mazE degrades faster than mazF, we built this constraint into the parameter sampling procedure. In addition, we considered the native mazEF system of E. coli. To prevent any potential interaction between the engineered and the native system, a mazEF deficient strain was created in the lab. This is, why we do not consider a native mazEF system in the model.

Click the button to see the detailed methods!

In Silico Kill Switch

We have built a mathematical model of pDusk and pDawn which describes the system behavior. Initial results from the lab, did not show a significant change in response for the pDusk system. However, the experimental results are in agreement with the computed data in the pDawn system. This leads us to suspect that either there was an irregularity with the dark conditions in the lab, or the model does not fully reflect the underlying biology of the system in the lower light intensity regions. We suggest that a parameter sensitivity and robustness analysis could be conducted to find the most likely cause of the mismatch to the data. We further conclude, that we achieved a good fit in Figure 12 to the data available in literature1.

We then studied the mazEF toxin-antitoxin integration to design our optogenetic kill switch. After analysis of this system, we found two parameter sets which satisfy the constraints. These parameters give us an idea on how to search the parameter space for further analysis to select functioning kinetic rates in collaboration with the wet-lab.
In addition, future work could show the effect of parameter sensitivity such as leaky promoters and their effect on the system to engineer the pDusk system, adjusting it to suit our needs. We also made the Matlab code available for future iGEM teams to build optogenetic tools.

Metabolic Modeling

Home Model Lab comic

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 m o l s u c r o s e L or above1. Supplementing an Apis mellifera (honey bee) colony with sugar water is a well established practice amongst beekeepers2. This is usually done with about 1 k g of table sugar for each liter of table sugarChemically speaking this is pure sucrose: C12H22O11. After heating and stirring, this ends up as a concentration of 625 g s u c r o s e L W a t e r . This concentration can also be defined as 1.82 m o l S u c r o s e L , 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?

Modeling 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 iJO13663 from the BIGG database. We modified the model to represent an E. coli in a sugar water-like environment using CobraPy4. We also changed the objective of the optimization such that the ATP-maintenance reaction should be maximized. Then we ran the model 1000 times whilst sequentially increasing the lower bound of the water efflux exchange reaction, the results of which can be seen in figure 17.

Figure 17: The relationship between the maximum ATP available for survival for an E. coli in a sugar-water environment and the incrementally increased water efflux.

The relationship between maximum ATP available for survival and water efflux is shown in Figure 17 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   m m o l g D W ⋅ h 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 17 we can start looking at the relationship between survival time and water efflux.

Figure 18: The relationship between the maximum ATP available for survival for an E. coli in a sugar-water environment and the theoretical survival time, given a constant water efflux over this time and a starting volume of 2.8e-13 grams5. The various coloured lines indicate water tolerance thresholds for the E. coli.

In Figure 18 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.

Beyond 90 Minutes

Figure 17 shows us that osmotic pressure alone can indeed have an effect on cell regulation and cell death, and from Figure 18 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. 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 previously published model beehave was adapted to include the effects 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 Disorder1,2. It consists of several modules which control aspects like foraging, Varroa mite dynamics, and colony growth as shown in Figure 19, and was extensively tested for robustness and realism1.

BeeT can be transported to the hive by adding it in sugar water or bee bread. Both these sources have their advantages and disadvantages. In this project we have used Escherichia coli as our model chassis, which can be applied to sugar water. Sugar water is supplemented to support a colony during honey harvest before winter, as a substitute for nectar3, 4. Supplementing honey bees with sugar water is a well established and familiar practice among beekeepers3. Alternatively, reconstructing the same system in Lactobacillus species would allow us to use BeeT in bee bread. Bee bread has some advantages over sugar water, as it it is more frequently transported to larvae. Bee bread is a combination of pollen, regurgitated nectar, and glandular secretions. To protect against harmful micro-organisms honey bees also inoculate bee bread with fermenting bacteria5. The fermentation process reduces the pH .

BeeT Module

Beehave is an open-source GNU licensed agent-based model utilizing NetLogo NetLogo is an open source, agent-based modeling platform widely used in education and research. , and consists of several interlocking modules which each model different aspects of the beehive. It models the wide variety of stresses affecting honey bees1. 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 colony dynamics, foraging, and a Varroa mite model as depicted in Figure 19. Two viruses are also included in the model: deformed wing virus (DMV), and acute paralysis virus (APV), both 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 honey and pollen stores. The Varroa mite mortality is determined by the amount of BeeT near larvae when Varroa mites emerge from brood cells. This in turn affects Varroa mite population levels in the hive, reducing virus loads in the hive, and allowing colony survival.

Figure 19: The honey bee colony model includes Varroa mite dynamics and virus epidemiology, agent-based foraging behavior with either pre-defined landscape definitions or a representation of local floral patterns. It is also possible to include weather patterns to more accurately model local conditions. Note that the model includes various interdependent mortalities and other parameters which are not included in this figure.

Model Considerations

Not everything is known about the BeeT model and as such, it was necessary to make certain assumptions. Some of these 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 sugar water, BeeT is exposed to the honey stomach of bees, honey stores, and brood cells, whereas bee bread encounters pollen and brood cells. 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 modeled this by using a saturating The chosen saturating function smoothly approaches one and is used to calculate mite mortality. function. We also assumed that the effect of BeeT on Varroa mite mortality is entirely determined by the amount of BeeT present at larvae when a brood cell is capped Nursing bees cover up brood cells to protect larvae, also called capping. . Finally, the beehave model is only able to model a single virus: DMV or APV. We are unable to model the combined effects of both viruses on a bee colony. For all analyses we used the DWV virus, as it is more harmful to honey bee colony survival than APV8.

Sugar Water: we know that, based on the experimental results obtained by metabolic modeling of E. coli, and it is able to survive in sugar water for at least 48 hours. Additionally, modeling predicts that if E. coli is capable of surviving 90 minutes, it can survive indefinitely. Additionally, we estimate that in sugar water there are approximately 10 6   c e l l s m l . This estimate is based on the assumption that we add saturated medium ( 10 8   c e l l s m l ) 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, since it would be too contaminated to consume. We based the hourly uptake rate of sugar water on a paper by Avni et al9, but it is possible that this is not the maximal uptake rate of sugar water. In the measured time period all of the sugar water was taken up by the colony, suggesting that the uptake rate is theoretically greater than the value proposed by Avni et al 9.

Bee Bread: we assume that the amount of Lactobacillus cells in bee bread is roughly equal to the number of C F U m l in yoghurt 10 . Artificial bee breads can be made using yoghurt and are readily accepted by bees10. Additionally, we assume that the removal rate of artificial bee bread is equal to 22.7 g d a y 4 . Based on Avni et al. we know that the uptake rate depends on the manner in which bee bread is applied9. We chose to base the uptake rate of bee bread on Brodschneider et al., since their experimental setup is similar to Avni et al4,10.

How Can Honeybees Survive?

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 nearby larvae. We are interested in how BeeT can best be applied given certain assumptions. As such we are interested in the best period to apply BeeT, either in autumn 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. Each is varied across a range of values. Additionally, for every combination we have five replicates to show the variance in our results due to stochastic processes. We divided the analyses into four treatments, each with the same three parameter ranges and five replicates per combination of parameters. This can be seen in Table 3.

Table 3: Structure of BeeT module analysis.
Period Treatment Range
1 April -
15 June

Spring 		Sugar water Varroa mite mortality, degradation in-hive and degradation outside-hive
Bee bread Varroa mite mortality, degradation in-hive and degradation outside-hive
1 Sep. -
15 Nov.

Autumn 		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 5169 parameter set per ranges. We compared the autumn period with the spring period for the different ranges and treatments. In each case the spring period was significantly more effective in combating Varroa mites than the autumn 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 winter. If the population is below 4000 worker bees the colony is presumed dead and the population is set to zero. To differentiate between the scenarios we looked at the mean population of worker bees after winter averaged over 10 years. Included in this average are the years where the bee population is zero. If the average is below 3000, we assume that the colony dies within a few years and consider this colony death. If the mean population after winter is between 3000 and 5000, the colony survives but its population is reduced due to Varroa mite pressure. Above 5000 worker bees, the colony is thriving and the initial infestation of Varroa mites rapidly dies off. 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 divide parameter sets in three categories: death, survival, and thriving. As can be seen in Table 4, the spring period is more effective for both sugar water and bee bread. Using bee bread is advisable when BeeT has a limited effect on mite mortality. If BeeT has a large impact on mite mortality it is more beneficial to use sugar water. This is because beekeepers are already comfortable applying sugar water to hives which reduces the risk of mistakes.

Table 4: Parameter sets are divided into three categories: colony death, survival, and thriving. If colonies can survive and thrive with higher degradation of BeeT (in-hive and outside the hive) and a lower effect of BeeT on Varroa mite mortality, it indicates a more effective treatment.
Period and treatment Colony death Colony survival Colony thriving
Sugar water, spring 6,6% 80,6% 12,8%
Bee bread, spring 0% 2,9% 97,1%
Sugar water, autumn 15,1% 80,7% 4,2%
Bee bread, autumn 0 57,6% 42,4%

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.

For sugar water, parameter sets were chosen which fall in the categories survival, and thriving. These are further analysed by examining the population dynamics of worker bees and Varroa mites over time. In Figure 20B and Figure 20C these are shown, respectively, for a colony barely surviving a Varro mite infestation, and a thriving colony, 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. For comparisons sake we also show in figure 20A with population dynamics when BeeT.

Figure 20: The honey bee population is shown in blue and the Varroa mite population in red. A: Colony rapidly declines when no BeeT is present. Starting population is 20 Varroa B: Colony barely survives Varroa mite infestation. Shows Varroa mite in red and worker bee population in blue. Starting population is 20 Varroa. C: Colony thrives regardless of Varroa mite infestation. Starting population is 20 Varroa mites and rapidly declines. D: Colony thrives regardless of heavy Varroa mite infestation. Starting population is 10,000 Varroa mites.

When the colony is thriving, the starting Varroa mite population rapidly declines and disappears as shown in Figure 19C. In the Figure you can see that the Varroa mite population is never able to establish itself. This indicates that the treatment is effective for a low population of initial Varroa mites. We further analysed a situation with a starting population of 10,000 Varroa mites which is shown in Figure 20D. 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.

How Was The Model Analysed?

The beehave model utilizes multiple number generators to calculate mortality amongst Varroa mites and bees. These 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.

Figure 21: Mean worker bee population over 10 years plotted against year round Varroa mite mortality at emergence from 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 compromize between simulation time Over 110,000 parameter sets were simulated. and reduction of variance.

We first separated the results of the analyses of the full BeeT module based on the three scenarios. The results can be seen in Figure 22. Within this 2D graph the axis are either degradation in the hive or degradation outside the hive, and impact of BeeT on Varroa mite mortality. Additionally, the mean population of worker bees surviving winter is color coded. We can see that the highest mean population is in the bottom right corner with the scaling factor affecting the results strongly. Additionally, degradation inside the hive has a greater impact on mean population than degradation outside the hive. This can be seen by comparing the top right corner of both panels.

Figure 22: 2D graph with either degradation in the hive or degradation outside the hive, and impact of BeeT on Varroa mite mortality on its axes. Results are color coded with mean population of bees after winter. The left panel has degradation outside the hive on its y-axis, degradation inside the hive is averaged out for each point. The right panel has degradation inside the hive on its y-axis, degradation outside the hive is averaged out for each point.

BeeT Beehave - Discussion and Conclusions

Currently it is unclear how effective BeeT will be at combating Varroa mites since, primarily, toxicity is difficult to measure in the lab and deploying BeeT in the field has many regulatory issues. For these reasons, we built 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: autumn and during spring. We chose autumn to analyse the effect of BeeT on the Varroa mite population since honey bee populations are at their most fragile during the winter following this period. This way we strengthen the honey bees since winter is coming. Additionally, it is normal for current beekeeper practices to supplement hives with sugar water during autmn and, thus, applying BeeT at this time 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 health6,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, supplementing sugar water at this time is already an accepted practice among 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 bees1. As such only a small fraction is transported to larvae. On the other hand, the protein contained in bee bread is crucial for the development of larvae and, therefore is mainly given to larvae to support their growth1. 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

We examined the effect of year round mite mortality by adding the following code to beehave:


which 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 the mite mortality ScalingVariable is large, the mortalityMites function approaches one.
The amount of Varroa mites emerging from a cell is calculated with the variable ‘healthyMitesInSingleCell’. Using a 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.

The initial code for 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 dependent 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 below, including how they are calculated, if applicable what their initial value is, units, in what procedure they are used, and any comments or references. The original variables of the beehave model can be found here.

Click the button to see the detailed methods!

References

    References for Population Dynamics

    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. Nazanin Saeidi, Mohamed Arshath, Matthew Wook Chang, Chueh Loo Poh, Characterization of a quorum sensing device for synthetic biology design: Experimental and modeling validation. Volume 103, 15 November 2013, Pages 91–99 ↩

    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. 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. 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. 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) ↩

    References for Optogenetic Kill Switch

    1. Ohlendorf, Robert, Roee R. Vidavski, Avigdor Eldar, Keith Moffat, and Andreas Möglich. 2012. “From Dusk till Dawn: One-Plasmid Systems for Light-Regulated Gene Expression." Journal of Molecular Biology 416(4):534–42. Retrieved (http://dx.doi.org/10.1016/j.jmb.2012.01.001) ↩

    2. Möglich, Andreas, Rebecca A. Ayers, and Keith Moffat. 2009. “Design and Signaling Mechanism of Light-Regulated Histidine Kinases." Journal of Molecular Biology 385(5):1433–44. Retrieved (http://dx.doi.org/10.1016/j.jmb.2008.12.017). ↩

    3. Klose, Cornelia et al. 2015. “Systematic Analysis of How Phytochrome B Dimerization Determines Its Specificity." Nature Plants 1(July):15090.↩

    4. Rausenberger, Julia et al. 2010. “An Integrative Model for Phytochrome B Mediated Photomorphogenesis: From Protein Dynamics to Physiology." PLoS ONE 5(5).↩

    5. Engelberg-Kulka, Hanna, Ronen Hazan, and Shahar Amitai. 2005. “mazEF: A Chromosomal Toxin-Antitoxin Module That Triggers Programmed Cell Death in Bacteria." J. Cell Sci. 118:4327–32.↩

    6. Limpert, Eckhard, Werner a. Stahel, and Markus Abbt. 2001. “Log-Normal Distributions across the Sciences: Keys and Clues." BioScience 51(5):341.↩

    7. Fasani, Rick a and Michael a Savageau. 2013. “Molecular Mechanisms of Multiple Toxin-Antitoxin Systems Are Coordinated to Govern the Persister Phenotype." Proceedings of the National Academy of Sciences 110(27):E2528–37.↩

    8. Mancinelli, A. L. (1994) The physiology of phytochrome action. In Photomorphogenesis in Plants, 2nd ed. (Edited by R. E.Kendrick and G. H. M.Kronenberg), pp. 211–269. Kluwer Academic Publishers, Dordrecht, The Netherlands.↩

    References for Metabolic Modeling

    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. 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. ↩

    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. Ebrahim, A., Lerman, J. A., Palsson, B. O., & Hyduke, D. R. (2013). COBRApy: constraints-based reconstruction and analysis for python. BMC systems biology, 7(1), 74. ↩

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

    References for Beehave

    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. ↩