A deeper understanding of system interactions can explain contradictory field results on pesticide impact on honey bees
The bee health modelThe conceptual model of the interactions of various stressors with honey bee health is described by the following system of ordinary differential equations (ODEs)$${{tau }_{{HB}}dot{x}}_{{HB}}= {-{delta }_{{HB}}x}_{{HB}}+{g}_{{TC}}left({x}_{{TC}}right)+{g}_{{VA}}left({x}_{{VA}}right)+{g}_{{VI}}left({x}_{{VI}}right) \ +{bar{f}}_{S,C}left({u}_{S},{u}_{C},{x}_{{TC}},{x}_{{VA}}right)+{bar{f}}_{P}left({u}_{P},{x}_{{TC}}right)+{underline{f}}_{{HB}}left({u}_{T}right)$$
(1)
$${{tau }_{{TC}}dot{x}}_{{TC}}={-{delta }_{{TC}}x}_{{TC}}+{g}_{{HB}}left({x}_{{HB}}right)$$
(2)
$${{tau }_{{VA}}dot{x}}_{{VA}}={-{delta }_{{VA}}x}_{{VA}}+{h}_{{VA}}left({{x}_{{HB}},x}_{{TC}},varepsilon {x}_{{VI}}right)+{underline{f}}_{{VA}}left({u}_{T}right)$$
(3)
$${{tau }_{{VI}}dot{x}}_{{VI}}={-{delta }_{{VI}}x}_{{VI}}+{h}_{{VI}}left({{x}_{{HB}},x}_{{TC}},{varepsilon x}_{{VI}}right)$$
(4)
for the state variables ({x}_{{HB}}) representing honey bee health, ({x}_{{TC}}) the stress due to toxic compounds (e.g., neonicotinoid insecticides), ({x}_{{VA}}) the stress due to parasites (e.g., V. destructor) and ({x}_{{VI}}) the stress due to pathogens (e.g., DWV). The system includes the effects of external inputs as sugar ({u}_{S}), pollen ({u}_{P}), absolute deviation from desired temperature ({u}_{T}) and sub-optimal temperature ({u}_{C}). All the inputs and possible parameters are non-negative; the coefficients (tau) denote the time constants; the coefficients (delta) denote the self-regulation parameters; (varepsilon) in the last two equations allows to account for pathogens that can ((varepsilon , > , 0)) or cannot ((varepsilon=0)) impair the immune system (through link m in Fig. 1). We assume that the functions (g) are smooth, bounded, positive, convex and decreasing to 0; the functions (bar{f}) are smooth, bounded, non-negative, concave and increasing with respect to (w.r.t.) (u) arguments (vanishing only when the first (u) argument vanishes) while convex and decreasing to 0 w.r.t. (x) arguments; the functions ({underline{f}}) are smooth, bounded, non-positive and decreasing (vanishing only when (u=0)); the functions (h) are smooth, bounded, positive, convex and decreasing to 0 w.r.t. the first argument while concave and increasing w.r.t. all the other arguments. For a detailed description of the various functions, together with a summary of the biological effects they account for and a reference to the conceptual model in Fig. 1, see Supplementary Table 3.Structural analysis of the bee health modelWe describe here the structural considerations and computations that yield the structural influence matrix for the honey bee health system.The structural influence matrix (M) is defined as follows. (M) is a symbolic matrix with entries ({M}_{{ij}}) chosen among: +,−,0,?, according to the criteria described below. Consider an equilibrium point (bar{x}) and a constant perturbation (u) applied on the (j)-th system variable (small enough not to compromise the stability of the equilibrium). The equilibrium value will be modified as (bar{x}+delta bar{x}). Consider the sign of the perturbation of the (i)-th variable, (delta bar{{x}_{i}}). Then ({M}_{{ij}}) = + if (delta bar{{x}_{i}}) always has the same sign as (u); ({M}_{{ij}}=) − if (delta bar{{x}_{i}}) always has the opposite sign as (u); ({M}_{{ij}}) = 0 if always (delta bar{{x}_{i}}=0); regardless of the system parameters. Conversely, if the sign does depend on the system parameters, we set ({M}_{{ij}}) = ?.In this section we prove that the influence matrix of the honey bee health system is structurally determined, i.e., there are no “?”‘ entries in (M).We start with the following proposition.
Proposition 1
Assume that a matrix
(J)
is Hurwitz stable (i.e., all its eigenvalues have negative real part) and has the sign pattern
$${sign}left(Jright)=left[begin{array}{cccc}- & – & – & -\ – & – & 0 & 0\ – &+& – &+\ – &+& 0 & -end{array}right]$$
Then, the sign pattern of
({adj}left(-Jright))
, the adjoint of
(-J)
, is
$${sign}left({adj}left(-Jright)right)=left[begin{array}{cccc}+& – & – & -\ – &+&+&+\ – &+&+&+\ – &+&+&+end{array}right]$$
Proof To prove the statement, we just change the sign of the first variable, hence we change sign to the first row and column of matrix (J). The resulting matrix (M) is such that$${sign}left(Mright)=left[begin{array}{cccc}- &+&+&+\+& – & 0 & 0\+&+& – &+\+&+& 0 & -end{array}right]$$We observe that (M) is a Metzler matrix, namely, all its off-diagonal entries are non-negative. Moreover, the matrix is Hurwitz stable. Then, we can proceed as in the proof of Proposition 4 in a previous report16. Given a Metzler matrix that is Hurwitz stable, its inverse has non-positive entries; hence, the inverse of (-M) has non-negative entries: ({left(-Mright)}^{-1}ge 0) elementwise. Moreover, we observe that(,M) is an irreducible matrix, i.e., there is no variable permutation that brings the matrix in a block (either upper or lower) triangular form. This implies that the inverse of (-M) has strictly positive entries: ({left(-Mright)}^{-1} , > , 0) elementwise. Also, stability implies that the determinant of (-M) is positive: ({det }left(-Mright) , > , 0). Then, ({adj}left(-Mright)={left(-Mright)}^{-1}{det }left(-Mright) > 0), hence the adjoint of (-M) is also positive elementwise. To consider again the original sign of the variables, we change sign to the first row and column of ({adj}left(-Mright)), and we get the signature above for ({adj}left(-Jright)).The next step is the characterization of the structural influence matrix, which corresponds to the sign pattern of the adjoint of the negative Jacobian matrix in Proposition 1.To this aim, we first consider the linearized system and write it in a matrix-vector form$$dot{x}left(tright)={Jx}left(tright)+{e}_{j}u$$where (dot{x}left(tright)) is the time derivative of the four-dimensional vector (xleft(tright)) and ({e}_{k}), (k={{{{mathrm{1,2,3,4}}}}}), is an input vector, constant in time, with a single non-zero component, the (k)-th, equal to 1, while the scalar (u , > , 0) is the magnitude of the input. We wish to assess the (i)-th component of (xleft(tright)), ({x}_{i}left(tright)={e}_{i}^{T}xleft(tright)). If (J) is Hurwitz, as assumed, the steady-state value of variable ({x}_{i}left(tright)) due to the input perturbation ({e}_{k}) applied to the equation of variable ({x}_{k}left(tright)) is achieved for$$0=Jbar{x}+{e}_{k}u,$$namely$${x}_{i}=-{e}_{i}^{T}{J}^{-1}{e}_{k}u,$$which implies that the sign of the steady-state value ({bar{x}}_{i}) of variable ({x}_{i}) due to a persistent positive input acting on the (k)-th equation has the same sign as ({(-{J}^{-1})}_{{ik}}), the (left(i,kright)) entry of matrix ({left(-Jright)}^{-1}). Since we assume Hurwitz stability, we have that ({det }left(-Jright)) is positive, hence the sign pattern of the inverse ({left(-Jright)}^{-1}) corresponds to the sign pattern of the adjoint, ({adj}left(-Jright)). In fact, ({adj}left(-Jright)={left(-Jright)}^{-1}{det }left(-Jright)).We next consider the nonlinear system under investigation, which we write in the form$$dot{x}left(tright)=fleft(xleft(tright)right)$$and without restriction we assume that the zero vector is an equilibrium point: (0=fleft(0right)). This condition can be always achieved, without loss of generality, by a translation of coordinates. We also consider a stable equilibrium: we assume that the linearized system at the equilibrium is asymptotically stable, namely its Jacobian (J), which has the sign pattern considered in Proposition 1 above, is Hurwitz. We also assume that a constant input perturbation of magnitude (u) is applied to the system, affecting the (k)-th equation, i.e.,$$dot{x}left(tright)=fleft(xleft(tright)right)+{e}_{k}u,$$and that the perturbation is small enough to keep the state in the domain of attraction of the considered equilibrium. Due to this perturbation, a new steady state (bar{x}left(uright)) is reached that satisfies the condition$$0=fleft(bar{x}left(uright)right)+{e}_{k}u$$To determine the sign of the new equilibrium components (bar{x}left(uright)), we consider this new equilibrium vector as a function of (u) in a small interval (left[0,{x}_{{MAX}}right]). Adopting the implicit function theorem yields$$frac{d}{{dx}}bar{x}left(uright)=-J{left(uright)}^{-1}{e}_{k}u,$$where we have denoted by (Jleft(uright)) the Jacobian matrix computed at the perturbed equilibrium (bar{x}left(uright)). Hence, for (u) small enough, the sign of the derivatives of the entries of the new, perturbed equilibrium are, structurally, the same as those in the (k)-th column of matrix (-{J}^{-1}). Since, by construction, (xleft(0right)=0), this is also the sign of the elements of vector (bar{x}left(uright)), for (u) in the interval (left[0,{x}_{{MAX}}right]).We have therefore proved that the original nonlinear system describing honey bee health admits the following structural influence matrix:$$left[begin{array}{cccc}+& – & – & -\ – &+&+&+\ – &+&+&+\ – &+&+&+end{array}right]$$System equilibriaThe results concerning the system equilibria were obtained through a standard analytical treatment of the nonlinear equations describing the equilibrium conditions of the system of differential Eqs. (1), (2), (3), (4). A detailed description of methods is reported in Supplementary Methods.Laboratory experiments using honey beesTo confirm the bistability of the system representing honey bee health as affected by multiple stressors, we used data from several survival experiments, carried out in a laboratory environment according to the same standardized method, over a 6-year period (Source data file).All experiments involved Apis mellifera worker bees, sampled at the larval stage or before eclosion, from the hives of the experimental apiary of the University of Udine (46°04′54.2″N, 13°12′34.2″E). Previous studies indicated that the local bee population consists of hybrids between A. mellifera ligustica and A.m. carnica62,63. Ethical approval was not required for this study.We considered experiments on the effect of the following stressors: infection with 1000 DWV genome copies administered through the diet before pupation, feeding with a 50 ppm nicotine in a sugar solution at the adult stage, exposition to a sub-optimal temperature of 32 °C at the adult stage. All experiments were replicated 3 to 13 times, using, in total, the number of bees reported in Table 1.For the artificial infection with DWV, we collected with soft forceps individual L4 larvae from the brood cells of several combs. Groups of 20–30 of such larvae were placed in Petri dishes with an artificial diet made of 50% royal jelly, 37% distilled water, 6% glucose, 6% fructose, and 1% yeast. 25 DWV copies per mg of diet were added or not to the diet according to the experimental group (note that a bee larva at this stage consumes about 40 mg of larval food per day, thus the viral infection per bee was 1000 viral copies). After 24 h larvae were transferred onto a piece of filter paper to remove the residues of the diet and then into a clean Petri dish, where they were maintained until eclosion. At the day of emergence, bees were transferred to plastic cages in a thermostatic cabinet, where they were kept until death. The DWV extract was prepared according to previously described protocols64 and quantified according to standard methods.For the treatment with nicotine, 10 µL of pure nicotine were added to 200 g of the sugar solution used for the feeding of the caged bees, to reach the concentration of 50 ppm.Finally, to expose bees to a 32 °C temperature, the plastic cages with the adult bees were kept in a thermostatic cabinet whose temperature was set accordingly.To monitor the survival of the adult bees treated as above, they were maintained from eclosion until death in plastic cages in a dark incubator at 34.5 °C (or 32 °C, according to the experiment), 75% R.H.; two syringes were used to supply a sugar solution made of 2.4 mol/L of glucose and fructose (61% and 31%, respectively) and water, respectively; dead bees were counted daily.All the results of these experiments are reported in Source data file.All experiments were carried out during the summer months, from June to September for 6 consecutive years. Previous data indicated that, in this region, virus prevalence increases along the active season starting from very low levels in spring and reaching 100% of virus-infected honey bees by the end of the summer; virus abundance in infected honey bees follows a similar trend28. For this reason, it can be assumed that bees sampled early in the season are either uninfected or they bear only a very low viral infection level, whereas bees sampled later in the season are likely to be virus-infected, bearing moderate to high viral infections. To confirm this assumption and identify a method for filtering our data according to viral infection, we assessed viral infection in a sample of bees from the untreated control group of each experiment, by means of qRT-PCR. According to standard practice, we assumed that Ct values below 30 are indicative of an effective viral infection, whereas Ct above that threshold are more likely in virus negative bees. As expected, we found that virus prevalence increases from June to September (Supplementary Figure 1a), in such a way that up to mid July only the minority of bees can be considered as viral infected (Supplementary Figure 1b). Therefore, we classified as “early” all the samples collected up to mid July and assumed that viral infection in those samples was low; on the other hand, samples collected from mid July till September were classified as “late” and we assumed that viral infection in those samples was high.qRT-PCR analysis of viral infection was carried out as follows. At the beginning of every experiment (i.e., at day 0), two to five bees for each replication were sampled in liquid nitrogen and transferred in a −80 °C refrigerator. After defrosting of samples in RNA later, the gut of each honey bee was eliminated to avoid the clogging of the mini spin column used after. The whole body of sampled bees was homogenized using a TissueLyser (Qiagen®, Germany). Total RNA was extracted from each bee according to the procedure provided with the RNeasy Plus mini kit (Qiagen®, Germany). The amount of RNA in each sample was quantified with a NanoDrop® spectrophotomer (ThermoFisher™, USA). cDNA was synthetized starting from 500 ng of RNA following the manufacturer specifications (PROMEGA, Italy). Additional negative control samples containing no RT enzyme were included. DWV presence was verified by qRT-PCR considering as positive all samples with a Ct value lower than 30. The following primers were adopted: DWV (F: GGTAAGCGATGGTTGTTTG, R: CCGTGAATATAGTGTGAGG65). 10 ng of cDNA from each sample were analyzed using SYBR®green dye (Ambion®) according to the manufacturer specifications, on a BioRad CFX96 Touch™ Real time PCR Detector. Primer efficiency was calculated according to the formula (E={10}^{left(-1/{{{{{{rm{slope}}}}}}}-1right)*100}). The following thermal cycling profiles were adopted: one cycle at 95 °C for 10 min, 40 cycles at 95 °C for 15 s and 60 °C for 1 min, and one cycle at 68 °C for 7 min.Individual survival and colony stabilityTo investigate how the death rate of forager bees affects colony growth, a compartment model of honey bee colony population dynamics was proposed50. This model showed that death rates over a critical threshold led to colony failure. Here we modified this model to include premature death of bees at younger age, as predicted by our model of individual bee health in the presence of an immuno-suppressive virus. We show that the critical threshold found in the previously published model50 becomes a decreasing function of the death rate of the younger individuals, so that premature death (and, in turn, immune-suppression) favors colony collapse.In more details, we first summarize the results of the previously published model50 where two populations (F) (forager) and (H) (hive) of bees are considered and where conditions are provided on the mortality (m) of (F) under which the whole population collapses: namely, mathematically stated, the system admits the zero equilibrium only. Here we extend the model partitioning (H) in two categories, (Y) (younger hive bees) and (O) (older hive bees), asintroducing an early mortality factor (n) for the young population, showing how such a factor worsens the collapsing condition.The previously published model50 concerns the interaction between hive bees (H) and forager bees (F) and is described by the ODEs$$dot{H}=Lfrac{H+F}{w+H+F}-Hleft(alpha -sigma frac{F}{H+F}right)$$$$dot{F}=Hleft(alpha -sigma frac{F}{H+F}right)-{mF}.$$Above, (L) is the queen’s eggs laying rate, (w) is the rate at which (L) is reached as the total population (H+F) gets large, (alpha) is the maximum rate at which hive bees become forager bees in the absence of the latter, (sigma) measures the reduction of recruitment of hive bees in the presence of forager bees and, finally, (m) is the death rate of forager bees (while the death rate of hive bees is assumed to be negligible).We first summarize the main results in terms of a threshold value for (m) in view of colony collapse, as our further analysis will follow a similar approach. All the parameters are assumed to be positive.The search for the equilibria of the above ODEs leads to the unique nontrivial equilibrium (beyond the trivial one)$$bar{H}=frac{L}{{mJ}}-frac{w}{1+J}$$$$bar{F}=Jbar{H}$$for$$J=Jleft(mright):=frac{alpha -sigma -m+sqrt{{left(alpha -sigma -mright)}^{2}+4malpha }}{2m}.$$Note that (J) is alway positive (and, moreover, it is independent of (L) and (w)). It follows that (bar{F}) and (bar{H}) have the same sign, so that the existence of the nontrivial equilibrium is equivalent to (bar{F}+bar{H} , > , 0). It is not difficult to recover that$$bar{F}+bar{H}=frac{w}{m}left(lfrac{1+J}{J}-mright)$$where (l:=L/w) is introduced for brevity. Then if (alpha le l) we get$$bar{F}+bar{H}=frac{w}{m}left(lfrac{1+J}{J}-mright)ge frac{w}{m}left(alpha frac{1+J}{J}-mright)=frac{w}{m}left(sigma+{mJ}right) , > , 0,$$with the last equality following from$$alpha -sigma frac{J}{1+J}-{mJ}=0,$$which in turn comes from annihilating the right-hand side of the second ODE and from using (J=bar{F}/bar{H}) while searching for equilibria. We conclude that, independently of (m), the colony never collapses if the recruitment rate (alpha) of forager bees is sufficiently low.Hence, we assume (alpha , > , l). Observe that$$bar{F}+bar{H}iff l , > , Jleft(m-lright)$$guarantees existence whenever (m) is sufficiently small, viz. (mle l). Assume then (m , > , l), so that the above condition reads$$J , < , frac{l}{m-l}$$leading to the threshold condition$$m , < , bar{m}:=frac{l}{2}frac{alpha+sigma+sqrt{{left(alpha -sigma right)}^{2}+4sigma l}}{alpha -l}$$by using the definition of (J), see Eq. (2) the previously published model50.A standard stability analysis shows that, assuming (alpha,m , > , l), the nontrivial equilibrium is (globally) asymptotically stable whenever it exists (positive), i.e., whenever (m , < , bar{m}). Otherwise, the only (globally) attracting equilibrium is the trivial one, corresponding to colony collapse (see Fig. 5 for the previously published model50 or Fig. 4 for (n=0)). In the mathematical jargon, the disappearance of the positive equilibrium, for (m) exceeding (bar{m}), is referred to as a transcritical bifurcation43.Now, in view of the outcome of the analysis of our model of individual bee health, we introduce a mortality term for the younger bees. As forager bees are recruited from adult hive bees, we divide the class of hive bees (H) in younger (Y) and older (O), assuming that the former die at a rate (n), while the death rate of the latter remains negligible according to the previously published model50. Obviously, (H=Y+O). The original ODEs are consequently modified as$$dot{Y}=Lfrac{H+F}{w+H+F}-Y$$$$dot{O}=left(1-nright)Y-Hleft(alpha -sigma frac{F}{H+F}right)$$$$dot{F}=Hleft(alpha -sigma frac{F}{H+F}right)-{mF}.$$Note that the sum of the first two equations above gives$$dot{H}=Lfrac{H+F}{w+H+F}-Hleft(alpha -sigma frac{F}{H+F}right)-{nY}.$$The new negative mortality term for younger hive bees, (-{nY}), models the fact that only the younger hive bees die prematurely while the rest of the dynamics is unchanged with respect to the original model.The search for equilibria soon gives$$bar{Y}=Lfrac{bar{H}+bar{F}}{w+bar{H}+bar{F}}$$from the first ODE above, so that the remaining two equilibrium conditions lead to$$bar{H}=frac{{L}_{n}}{{mJ}}-frac{w}{1+J}$$$$bar{F}=Jbar{H}$$for the same (J) originally defined and ({L}_{n}:=Lleft(1-nright)) (note that (nin left({{{{mathrm{0,1}}}}}right)), and the case (n=0) brings us back to the original model). From this point on the analysis is the same as that previously summarized for the original model, but for replacing (L) with ({L}_{n}) and (l) with (l:=lleft(1-nright)). Consequently, by assuming (alpha,m , > , {l}_{n}) (which is less restrictive when (n , > , 0)), the threshold condition (m < bar{m}) becomes$$m , < , bar{m}left(nright):=frac{{l}_{n}}{2}frac{alpha+sigma+sqrt{{left(alpha -sigma right)}^{2}+4sigma {l}_{n}}}{alpha -{l}_{n}},$$which clearly returns the original threshold condition when (n=0). Since$$frac{dbar{m}}{{dn}}left(nright) , < , 0$$as it can be immediately verified, it follows that the critical value for (m), (bar{m}left(nright)), beyond which the colony system admits only the zero equilibrium, i.e., the transcritical bifurcation value, decreases with (n) (Fig. 4). We thus conclude that colony collapse is favored by the premature death of younger hive bees, possibly caused by a virus impairing the immune system as shown by the analysis of our model of individual bee health.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
