Model | UNILausanne

Motivation

Our project focuses on limiting the spread of the invasive Asian hornet, which poses a significant threat to ecosystems across Europe. Our solution is addressed primarily to beekeepers, whose beehives are under constant attack by these aggresive species. To address this issue we have created a trap to deliver an RNA interference-based living insecticide specifically targeting the Asian hornet. This trap contains a genetically engineered bacterium (Lactococcus lactis), gut commensal of the Asian hornet, producing short hairpin RNA (shRNA). These shRNA molecules silences an essential gene (chitin synthase) of the Asian hornet larvae thereby disrupting its development.
Go to Description Page

During the course of our project, we actively questioned our initial assumptions and asked ourselves whether targeting larvae or adult workers would be the most efficient at fighting Asian hornets and whether our dsRNA should be delivered directly in the form of an shRNA or through bacteria producing the shRNA, that could colonise the gut of the Asian hornet larvae and spread among the colony thereby controlling their population.

Since answering both questions would require real-life experiments, which take a lot of time and expertise that are not directly available to us, and animal testing (which we want to limit), we decided to create a virtual simulation of a bee population under pressure from the Asian hornet invasion. Together with our assistant Dr. Gábor Hollo, an expert in modelling, we built upon an existing mathematical model (Romero-Leiton et al. 2022) of the bee population, to which we added new terms:

Hornet predation rate on the bee population
Hornet population, modelled after the bee population model, to which we added:

an infection rate by our engineered bacteria
a cross-infection rate from mature to immature hornets

On this page, you can find out in more detail how we designed our model and what results it gave us.

Jump to design section Jump to results section

Model Description

To create our model, we chose to extend an existing model on bee population (Romero-Leiton et al. 2022), which is based on two populations: one describing the mature bee population (B), and the other describing the immature stage (b).

We extended the bee population model with three new populations to simulate interactions with the Asian hornet population. The three new populations describe:

Mature hornets (H)
Immature hornets (h)
Hornets infected by our engineered bacteria (I)

This extended model enables us to capture the dynamics of the bee population over time, accounting for the pressures exerted by the Asian hornet.

See Table 1 for more details.

Model Design

To design the bee-hornet model, we followed a 6-step modeling cycle based on an introductory class on mathematical models given by Sara Mitri, a professor and researcher at the University of Lausanne, who studies microbial communities (The Mitri Lab – Evolutionary Ecology of Microbial Communities, s. d.).

This cycle consists of the following steps:

Design

Why do we need a model ?
What elements do we need in our model ?
What does our model look like ?
What are our equations ?

Analysis

How did we analyse our model?
What makes our model reliable? Did it help us answer our questions?

Results

How did we interpret our results?

1. Why do we need a model?

We needed a model to help us determine which strategy would be more efficient: targeting larvae or adult workers. Additionally, the model would predict the minimal infection rate at which our engineered bacteria should spread within the Asian hornet population to reduce its pressure on the honeybee population.

In short, we aimed to create a model that would provide insights such as:

“It is better to target larvae rather than adults!” or the reverse.
“Delivering dsRNA through a gut-commensal bacterium is better than directly feeding dsRNA to hornets!” or the reverse.

2. What elements do we need in our model?

The main elements needed in our model include:

A bee population,
A hornet population,
Mature and immature hornets,
An infected population of hornets,
An infection rate by our engineered bacteria,
A cross-infection rate between mature and immature hornets,
Factors of population decrease, such as death rates due to natural causes, stress, predation by Asian hornets, and infection,
Factors of population increase, such as production rates of the immature stage and maturation rates.

These points can be divided into two categories:

Populations that we want to simulate to observe their change over time, implemented as variables depending on time (Table 1).
Factors we will use to simulate an increase or decrease in population size, implemented as parameters (Table 2).

3. What does our model look like?

Do you remember what physics teachers would usually say to trick you into believing that you could easily solve their problems? “Start by making a drawing of all the information you have before solving any equation.”

Well, we did that and it worked!

We represented our populations (variables) with circles, and connected them with incoming and outgoing arrows symbolizing increases and decreases in population.

4. What are our equations?

To quantitatively describe bee and hornet populations’ change over time, we wrote five ordinary differential equations (ODE), one for each variable:

After considering all the above design aspects, we finally created our model:

In summary, our model builds upon an existing bee population model (Romero-Leiton et al. 2022) by incorporating three additional populations to represent the interaction with Asian hornets. These new populations (mature hornets (H), immature hornets (h), and hornets infected by our engineered bacteria (I)) allow us to better simulate the impact of hornets on bee populations over time.

If you are interested in our model in more detail, check out how we coded our model.

GitLab

Fun Fact: Even though we're already drowning in a sea of programming languages—HTML, CSS, and Java for our wiki; R, Bash, and Python for software development—we thought, "Why not add a little MATLAB to our chaos?" Because who doesn’t love getting their hands even messier?

Analysis

How did we analyze our model?

To analyze our model, we used MATLAB to solve equations and run simulations (GitLab link).

We first chose to independently solve the 2-dimensional model of bee population (eq. 1-2), since equations 3, 4, and 5 are not dependent on any bee variable.

We solved this system to find H*, the critical hornet concentration at which the bee population collapses. At this critical hornet concentration, the fixed point corresponding to a zero bee population becomes stable. After determining the fixed points, we performed a linear stability analysis to check their stability.

But why?

Answer: In a one dimensional representation, the derivative can be represented as the slope of a curve. When the time derivative reaches zero, the result of the function does not change anymore and the system has reached a stationary state. We follow the same principle to assess the stability of the fixed points in higher dimensions.

          
          
          %% We started by initializing our parameters that we mentioned in the above table for the model:
          %% Parameters
          p.beta = 1;     % Immature reproduction rate in [0.1,1]
          p.omega = 1;    % Adult maturation rate in [0.91,1]
          p.mub = 0.3;    % Natural immature death rate in [0.1,1]
          p.muB = 0.3;    % Natural mature death rate in [0.29,1]
          p.sigma = 0.1;  % Adult bee death rate from a stressful factor [0.07,0.82]
          
          p.k = 0.1;      % infection rate
          p.kh = 100;     % adult-immature infection rate
          p.delta = 1;    % infected death rate
          p.kB = 1;       % hornet eating bees
          
          %% SECTION 1
          %% bee subsystem
          %% POPULATION CHANGE OVER TIME BASED ON REPRODUCTION, MATURATION AND DEATH.
          % symbolic variables
          syms beta omega mub muB sigma kB b B H
          
          % the rhs of the ODE
          f = [beta*B/(1+B)-omega*b-mub*b;
              omega*b-muB*B-sigma*B-kB*H*B];
          
          assume([beta omega mub muB sigma kB] > 0)
          
          % solve the equation
          sol = solve(f==0,[b,B]);
          
          % Substitutes the parameters into the solution while setting the hornet population H to zero
          % Analyzing the bee subsystem without hornets to identify their intrinsic population dynamics without predation pressure
          sol.b = subs(sol.b,[beta omega mub muB sigma kB H],[p.beta p.omega p.mub p.muB p.sigma p.kB 0]); 
          sol.B = subs(sol.B,[beta omega mub muB sigma kB H],[p.beta p.omega p.mub p.muB p.sigma p.kB 0]);
          
          %% Reaction Kinetics: Simulating population dynamics of bees over time
          % Time integration limits
          tmin = 0;
          tmax = 1e2;
          
          % Initial conditions from a hornet-free fixed point
          c0 = [double(sol.b(1)) double(sol.B(1)) 1e-3 0 0];
          
          % solve the equation
          options = odeset('RelTol',1e-8,'AbsTol',1e-10);
          [t,c] = ode15s(@(t,c) reactions_bee_hornet(t,c,p),[tmin,tmax],c0,options);
          
          % plot the results
          hf = figure;
          plot(t,c,'LineWidth',2)
          legend('Immature bee','Mature bee','Immature hornet','Mature hornet','Infected hornet','Location','northeast')
          
          % labels
          xlabel('\\bf\\boldmath$t$ / week','interpreter','latex','Fontsize',18)
          ylabel('\\bf Species population','interpreter','latex','Fontsize',18)
          
          % settings for the given axis
          set_ax(gca);
          
          % save figures
          exportgraphics(hf,['trajectory_kh_' int2str(p.kh)  '.pdf'],'Resolution',300)

To solve the two-dimensional ODE system of the bee population, we derive equations 1 and 2 to find the Jacobian matrix (J) (eq. 6) of the system:

In equation 6, H* represents a stationary hornet population.

When the bee population collapses, the system should reach a stable fixed point at (b, B) = (0, 0), when there are neither mature nor immature bees left.

We did not calculate the exact eigenvalues of the Jacobian matrix (eq. 6), as our focus was solely to determine the stability of the fixed points. To speed up the simulations, we ensured the stability of the (0, 0) fixed point with the Routh-Hurwitz criterion:

The stability criterion for the fixed points is that the real parts of the eigenvalues must be negative. For a two-dimensional system, it can be fulfilled if the following two criteria are satisfied:

The trace of the Jacobian (eq. 7) must be negative,
The determinant of the Jacobian (eq. 8) must be positive.

Since all parameters used in the model are positive (see Table 2), equation 7 is always satisfied.

From equation 8, we can calculate at what concentration the hornet population preys too much on bees and induces the collapse of the bee population (see eq. 9):

Based on our selected parameters (see Table 2), the bee population collapses when 36.92% of the population of Asian hornets prey on bees.

To further analyze our model, we examined it numerically using Matlab (see GitLab link), since it cannot be solved analytically.

We varied the following two parameters to assess their effect on the model (see Fig. 2):

Cross-infection rate (k_h), varied between 0 (no transmission from infected adult workers to larvae) and 100 (population time)^-1,
Infection rate (k), varied between 0 (no infection by our engineered bacteria) and 1 (time)^-1.

      
        
      %% This section defines the full system of ODEs that includes both bees and hornets, 
      %% capturing more complex interactions such as infections and predation. 
      %% The system now includes immature and mature hornets as well.
      
      %% The ODE system
      syms beta omega mub muB sigma k kh delta kB b B h H I
      
      assume([beta omega mub muB sigma delta kB] > 0)
      assume([k kh] >= 0)
      % assume([b B h H I] >= 0)
      
      % The rhs of the ODE
      f = [beta*B/(1+B) - omega*b - mub*b;
          omega*b - muB*B - sigma*B - kB*H*B;
          beta*H/(1+H) - omega*h - mub*h - kh*I*h;
          omega*h - muB*H - sigma*H - k*H;
          k*H - delta*I];
      
      % Substitute the given parameters
      %% In this part, the model calculates equilibrium points (fixed points) for both bee 
      %% and hornet populations, representing stable states where population sizes remain 
      %% constant over time. Understanding these fixed points helps to identify conditions 
      %% under which populations can coexist, thrive, or face extinction due to predation, 
      %% infection, or environmental stress.
      
      f = subs(f,[beta omega mub muB sigma delta kB kh],[p.beta p.omega p.mub p.muB p.sigma p.delta p.kB p.kh]);
      
      %% Calculating the fix points
      
      % Solve the equation
      sol = solve(f==0,[b,B,h,H,I]);
      
      % Number of fixed points
      nfix = numel(sol.b);

Figures 3A, 3B, 3C show a scenario where larvae are not infected (not targeted by our RNA interference), because transmission rate from infected adult workers to larvae is set to 0 (k_h = 0 (population time)^-1). Figures 3D, 3E, 3F show a scenario where larvae are infected (k_h = 100 (population time)^-1).

Figures 3A, 3D: Mature bee population as a function of infection rate (k); Figures 3B, 3E: Mature hornet population as a function of infection rate (k). The solid blue line represents the stable fixed point, the dashed red lines represent the unstable fixed point. The dotted black line represents the maximal tolerable hornet concentration (0.3692), since bee population collapses at this concentration.

In Figures 3C, 3F, example trajectories of each sub-population are plotted over time, starting without mature hornets and with 0.1% of immature hornets, simulating the beginning of an Asian hornet invasion. These trajectories were calculated with a 0.1 (week)^-1 infection rate, under which the bee population collapses if there is no cross-infection in the hornet population.

What makes our model reliable? Did it help us answer our questions?

Even though a model is a simplification of reality and may not be accurate in specific details, it should still make biological sense in order to be reliable. We ensured that our model behaves realistically by solving it for the stable fixed point (0,0), which means that when neither immature nor mature bees are present in the population, the only possible direction of the system is to stay at (0,0).

We chose to stick to a relatively “simple” model, as we believe its complexity is sufficient to answer our questions. However, if we had more time, we would still improve it by confirming the behaviour and dynamics of Asian hornet populations with an expert, since there is currently limited knowledge on this invasive species.

From Figure 3, we can find information to answer our 2 initial questions:

Plotting bee and hornet populations as a function of the infection rate (k) indicates that delivering shRNA through our engineered bacteria and gut colonisation helps keep the critical hornet concentration below the bifurcation value, represented by a dotted black line at 0.3692.
To understand if targeting larvae or adult workers would be more efficient, we compared 2 scenarios:
- one where there is no cross-infection of larvae by infected adults (Figures 3A, 3B, 3C),
- the other where there is cross-infection (Figures 3D, 3E, 3F).

In Figures 3B and 3E, the hornet population drops under the critical hornet concentration (H*) threshold as the infection rate (k) of adult workers increases. Without cross-infection (Fig. 3B), the hornet population reaches the threshold for an infection rate of 0.19, while with cross-infection (Fig. 3E), it reaches the threshold for a much lower infection rate, k = 0.01, ten times smaller than without cross-infection.

This suggests that we would not need a high number of hornets to be infected initially by our engineered bacteria if they spread the infection to the larvae at a sufficiently high rate.

Results

In Figures 4B and 4E, the hornet population drops under the critical hornet concentration (H* = 0.3692) threshold as the infection rate (k) of adult workers increases. Without cross-infection (Fig. 4B), the hornet population reaches the threshold for an infection rate of 0.19, while with cross-infection (Fig. 4E), it reaches the threshold for a much lower infection rate, k = 0.01, which is ten times smaller than without cross-infection.

Compared to Fig. 4A (no cross infection), cross-infection (k_h) pushes the bee population upwards (Fig. 4D), towards survival, and reduces the infection rate needed to avoid extinction (bee population at 0).

Compared to Fig. 4B (no cross infection), cross-infection (k_h) pushes the hornet population downwards (Fig. 4E), and reduces the infection rate needed to limit the hornet population under the critical hornet concentration threshold (dotted black line, H* = 0.3692).

This suggests that we would not need a high number of hornets to be infected initially by our engineered bacteria if they spread the infection to the larvae at a sufficiently high rate. Based on these results, we were convinced of the quality of our approach: we found that, theoretically, a solution exists that would allow us to reduce the pressure on beehives by infecting a relatively small number of adult hornets. This is in line with our attempt to target the Asian hornet while minimizing the harm on local fauna, and justifies the various specificity layers of our strategy. The model highlights that a crucial factor for the efficient eradication of the invasive species is the cross-infection rate. We tried to maximize this rate by engineering a commensal of the Asian hornet, in the hope that the bacteria proliferates and spreads in the hornet nest.

This evidence supports our idea that, within the context of our project, targeting larvae instead of adults and delivering shRNA through our engineered bacteria rather than administering shRNA directly are the most efficient approaches to limit the spread of the Asian hornet.

Over the course of our project, many experts have had a positive perspective on our idea of using an engineered bacterium to infect the Asian hornet, although the use of GMOs (genetically modified organisms) is currently highly controversial in many countries. This has led us to value mathematical modeling as a significant tool to evaluate efficiency and risks of GMO-based pest control strategies, providing more insights to the debate.

For instance, models could illustrate the consequences of failing to control the spread of an invasive species, by quantifying the speed and scale of damage caused by these species to agriculture and local biodiversity. Modeling could estimate the minimal amount of GMO required for effective biocontrol and predict how quickly invasive species might be eradicated thanks to GMO-based methods. Modeling could also compare the potential environmental impact of GMOs with the damage caused by invasive species if they are allowed to spread unchecked. In pest-management fields, mathematical models could provide guidance to decide on whether GMO use would be less harmful than ignoring invasive species.

References

Romero-Leiton, Jhoana P., Alejandro Gutierrez, Ivan Felipe Benavides, Oscar E. Molina, and Alejandra Pulgarín. 2022. “An Approach to the Modeling of Honey Bee Colonies.” Web Ecology 22 (1): 7–19. https://doi.org/10.5194/we-22-7-2022.
Russell, Stephen, Andrew B. Barron, and David Harris. 2013. “Dynamic Modelling of Honey Bee (Apis mellifera) Colony Growth and Failure.” Ecological Modelling 265 (September): 158–69. https://doi.org/10.1016/j.ecolmodel.2013.06.005.
“The Mitri Lab – Evolutionary Ecology of Microbial Communities.” Accessed September 30, 2024. https://wp.unil.ch/mitrilab/.
Torres, David J., Ulises M. Ricoy, and Shanae Roybal. 2015. “Modeling Honey Bee Populations.” PLOS ONE 10 (7): e0130966. https://doi.org/10.1371/journal.pone.0130966.
Winfree, Rachael. 2010. “The Conservation and Restoration of Wild Bees.” Annals of the New York Academy of Sciences 1195 (1): 169–97. https://doi.org/10.1111/j.1749-6632.2010.05449.x.