A Mathematical Model for Biocontrol of the Invasive Weed Fallopia japonica.

We propose a mathematical model for biocontrol of the invasive weed Fallopia japonica using one of its co-evolved natural enemies, the Japanese sap-sucking psyllid Aphalara itadori. This insect sucks the sap from the stems of the plant thereby weakening it. Its diet is highly specific to F. japonica. We consider a single isolated knotweed stand, the plant's size being described by time-dependent variables for total stem and rhizome biomass. It is the larvae of A. itadori that damage the plant most, so the insect population is described in terms of variables for the numbers of larvae and adults, using a stage-structured modelling approach. The dynamics of the model depends mainly on a parameter h, which measures how long it takes for an insect to handle (digest) one unit of F. japonica stem biomass. If h is too large, then the model does not have a positive equilibrium and the plant biomass and insect numbers both grow together without bound, though at a lower rate than if the insects were absent. If h is sufficiently small, then the model possesses a positive equilibrium which appears to be locally stable. The results based on our model imply that satisfactory long-term control of the knotweed F. japonica using the insect A. itadori is only possible if the insect is able to consume and digest knotweed biomass sufficiently quickly; if it cannot, then the insect can only slow down the growth which is still unbounded.

Introduction

The Japanese knotweed Fallopia japonica has been present in the UK since 1825. It was originally introduced by the Victorians as an ornamental plant, but soon escaped from their gardens. The plant can grow extremely quickly, forming very dense thickets which can shade out native plants and offer a poor habitat for native insects, birds and mammals, can grow in any soil type, and can penetrate concrete and cause damage to paved areas. DEFRA estimates the total cost of control of the weed in the UK to be £1.56 billion (see DEFRA 2003). One of the main problems with the plant is that it has no natural predators in the UK and therefore does not compete fairly with native species. Swansea is one of the worst affected areas, and it has been estimated that the current infestation there would cost 9.5 million to treat completely (see Cardiff 2006). The plant is also found in most of the states of the USA and several Canadian provinces.

The main quantifiable cost is herbicide treatment, but the weed has an extensive root system and can penetrate paving, tarmac and asphalt. It can damage flood defence structures and archaeological sites, and its removal can significantly add to the costs of preparing development sites. The weed is very difficult to eradicate and efforts at removal involving disturbing the soil can make matters worse. It forms very dense thickets that reduce native species diversity and are of little value to wildlife. It spreads easily: fragments of stem and rhizome are easily conveyed along rivers to distant sites, and the fibrous stems are slow to decompose. Reproduction of Japanese knotweed in the UK is currently by vegetative means only (regeneration from root material), and dispersal is due primarily to rhizome fragments being washed downstream in rivers, or the transportation by humans of soil containing rhizome fragments. Sexual reproduction does not occur in the UK as only the female plant is present there. In the native habitat, insect pollination, sexual reproduction, and wind dispersal of seed can also be contributing factors. Experimental work in the USA by Forman and Kesseli (2003) suggests considerable implications for management of the weed, should sexual reproduction start to become important.

Currently, control of F. japonica is primarily by chemical means using a herbicide such as glyphosate. This carries risks of contamination of rivers and streams where many knotweed areas are found, and in any case knotweed is in such rapid expansion that increased chemical use is not a long-term option. The shoots are edible by sheep, goats, cattle and horses and therefore satisfactory control is already achievable in grazing locations.

In this paper, we propose a mathematical model for the biocontrol of F. japonica using one of its co-evolved natural enemies, the Japanese sap-sucking psyllid Aphalara itadori. This insect sucks the sap from the stems of the plant thereby weakening it. After extensive experimental trials by the agricultural research organisation CABI (2007), in 2010 DEFRA in the UK approved the release of this insect as a biocontrol agent for F. japonica, initially at a handful of closely monitored sites. The insect was selected, among numerous species of plant-eating insects and fungi, because its diet is highly specific to F. japonica. Biological control has the advantage of being environmentally sound and more sustainable than extensive herbicide use, but the strategy will not result in complete eradication. The most that can be anticipated is that an equilibrium is reached in which the knotweed is still present but at an acceptable level. Such a strategy has been highly successful in controlling the Salvinia weed in Sri Lanka which was reduced by 80 % in four years by the introduction of the weevil Cyrtobagous salviniae, and in Kerala, India (Jayanth 1987). In the UK it seems likely that all knotweed is of one single clone (Hollingsworth and Bailey 2000), and the absence of sexual reproduction makes the plant a good candidate for biological control.

The model we develop in this paper is for a single isolated knotweed stand, and the plant’s size is described by time-dependent variables for total stem and rhizome biomass. Spatial effects are not modelled explicitly. However, the stems only grow to a certain height (about 3 m) and the roots only to a certain depth, also about 3 m. Therefore, unlimited growth of plant biomass, one of the predictions of our model under certain conditions, can only mean spatial spread of the stand. As far as the insects are concerned, it is the larvae of A. itadori that do the most damage to the plant and so the insect population is broken down into numbers of larvae and adults, using a standard stage-structured modelling approach. We assume that only larvae eat the stems (more precisely, the sap from the stems). Information on A. itadori can be found in Shaw et al. (2009). Eggs are laid on the upper surfaces of leaves, and the egg to adult duration is around 33 days (knotweed can grow one metre in this time). Insects can overwinter on sheltering plants. Neither larvae nor adults of A. itadori will eat the roots. Another species tested by CABI, the endoclyta moth, does eat the roots but was rejected as a biocontrol agent because it was found to also attack another unrelated plant (CABI 2007).

In this paper we want to give due consideration to the fact that the knotweed plant species exhibits unusual growth behaviour due to the fact that it is an invasive species that grows without the usual inhibiting factors of predation, or even lack of physical space. Given enough land, there is no limit to the extent to which it can spread over the time scale of interest and therefore even exponential growth is not unrealistic. The reader will notice that the way we model stem and rhizome growth does not provide for a carrying capacity and in fact the production rates for these quantities are linear. This is intentional and, far from being a simplification, it actually makes the model remarkably resistant to analysis that would usually lead in a straightforward way to properties such as boundedness of solutions, extinction or persistence using well-developed theories. The overall model is nonlinear, due to the fact that although the number of knotweed stems may grow without bound, each stem only grows to a finite height and therefore has a finite carrying capacity for eggs. This is modelled using a nonlinear function that represents the egg laying rate per stem, which levels off and may even decrease as the number of adults per stem increases due to crowding effects.

It turns out that the dynamics of the model depends mainly on a parameter h in our model, which measures how long it takes for an insect to handle (digest) one unit of F. japonica stem biomass. If h is too large then the model does not have a positive equilibrium and the plant biomass and insect numbers both grow together without bound, though at a lower rate than if the insects were absent. On the other hand, if h is sufficiently small, then the model possesses a positive equilibrium which appears to be locally stable. The structure of the model and its lack of monotonicity properties make it remarkably resistant to analytic study and therefore the analysis is based mostly on linear theory and backed up with numerical simulation. The overall take home message is simple: satisfactory long-term control of the knotweed F. japonica using the insect A. itadori is only possible if the insect is able to consume and digest knotweed biomass sufficiently quickly; if it cannot, then the insect can only slow down the growth which is still unbounded.

Model Derivation

For a complete list of notation, and the values given to parameters later in this paper, see Table 1. The state variables are A(t), L(t), S(t) and R(t), where A(t) and L(t) are the numbers of adult and larval A. itadori on a particular knotweed stand at time t, S(t) is the total stem biomass of the knotweed stand, and R(t) is total rhizome biomass. We begin by deriving equations that govern the evolution of these variables, starting with the larval A. itadori. Let a denote age and l(t, a) denote the age density of the larvae. Let (the subscript p standing for psyllid) denote the developmental time of A. itadori from egg to adult. Using a standard age-structured modelling approach, we may writewhich simply expresses the fact that larval A. itadori are lost at a per-capita rate that depends on the total stem biomass. The function will be decreasing, since larval mortality will be greater if there is a shortage of stems, and hence of sap, and is assumed to approach a positive limit as its argument tends to infinity, that limit being the natural per-capita mortality when there is no shortage of sap. We also write .

Table 1

Parameters used in the model and simulations

Symbol	Meaning	Value	Reference
τ_p	Development time of A. itadori from egg to adult (maturation time)	32.2 ± 0.5 days (28–42 days)	Shaw et al. (2009)
τ_s	Life-span of a stem	150 days	Approximation
k_l	Additional per-capita larval A. itadori mortality when sap is in short supply	k_l=1-m
m=μl∞	Natural per-capita larval A. itadori mortality	0.0205 per day	Shaw et al. (2009)
	When sap is plentiful
P	Maximum number of eggs per unit time that an individual adult can lay	23.9798 eggs per adult per day	Shaw et al. (2009)
q	Crowdedness constant	1	Approximation
μ_a	Per-capita mortality of adult A. itadori	0.0806 per day	Shaw et al. (2009)
μ_s	Per-capita natural mortality of stem biomass	0.01 per day per unit	Approximation
e	Resource (sap) encounter rate	1	Approximation
h (h_crit)	Handling (digestion) time per unit biomass consumed	The number of days per unit per larva	(To be estimated)
σ	Fraction of encountered food biomass that a larvae ingests	1	Approximation
k_s	Birth rate for stems	0.1 unit per day	Approximation
k_r	Birth rate for rhizome	0.1 unit per day	Approximation
μ_r	Per-capita loss rate of rhizome biomass	0.01 per day per unit	Approximation
A^0	Initial number of adult A. itadori	50	Approximation
R^0	Initial quantity of rhizome biomass	500 units	Approximation

For define . Differentiation shows thatThereforeand, setting ,Now, l(t, 0) is the birth rate of A. itadori, which should depend on the population of adult A. itadori and on the biomass of the knotweed stand available. Note that an individual stem can only accommodate a certain number of eggs so that intra-specific competitive effects apply at the level of the individual stem. Due recognition of this fact is crucial, since an important aspect we wish to study is the possibility that the knotweed stand as a whole can grow without bound, and therefore so can the number of larval and adult A. itadori. To reflect such an important competitive effect, we let be the egg laying rate per stem which is a function of the number of adult A. itadori per stem. We will keep as general as possible, but it is useful to keep in mind the frequently used prototype where (the well-known Nicholson’s blowflies birth rate), and where P is the maximum number of eggs per unit time that an individual adult can lay. This choice models the idea that the egg laying rate per stem increases nearly linearly with the number of adults per stem if the latter is fairly low, but reaches a maximum and then rapidly drops off if the number of adults per stem is very high, since then the issue of available space on the leaves for oviposition becomes important. The birth rate of A. itadori is therefore given byIf is chosen as , the implication is that the overall egg laying rate (i.e. for the whole knotweed stand) is . This gives the expected formula PA(t) if there are very few adults per stem, and a much lower value if A(t) / S(t) is very large. Note also that if which again is correct because if there are no stems there are no leaves and no oviposition sites, and the factor measures the impact of crowdedness on egg laying.

Using (3), (2) becomesand the total number of larvae isDifferentiating this yieldsThe last term in (6) is the rate at which larvae mature into adults. We may therefore immediately write down a differential equation for the number A(t) of adult A. itadori:where is the per-capita mortality of adult A. itadori.

We now derive an equation for the biomass S(t) of a stand. Shoots appear from the ground and rapidly form canes which are fully grown by the Summer and die back after that. The stems have their sap consumed by the larvae of A. itadori, and we assume this is the most important factor involved in stem damage, but we also allow for natural mortality of stem biomass, at a per-capita rate . Letting s(t, a) be the age density of stem biomass, we may suppose that stem loss is governed bywhere the first term in the right-hand side represents loss of stem biomass due to predation by A. itadori larvae. This is modelled essentially in terms of a Holling type II functional response. If stem biomass is very high, the sap consumption rate per larva would have to reach a plateau [the reader may wish to look ahead to see how the functional response eventually appears in the equation for S(t); see (12)]. In (8), the Holling type II functional response has been implemented with commonly used notation where e is the resource (sap) encounter rate, h is the handling (digestion) time per unit biomass consumed and is the fraction of encountered food biomass that the larva ingests. Using similar calculations to those that led to the expression for l(t, a) in (4), we can show thatNow, s(t, 0) is the rate of production of stem biomass, i.e. of new shoots. These originate from the roots and so we let be the production rate of stem biomass, where the subscript s stands for stems, and writeUsing (10) in (9), we find that the total stem biomass is, after a change in the variable of integration, given bywhere is the life-span of a stem, i.e. the time between the appearance of the new shoot and the time when the stem dies back. Differentiating (11) shows thatThe first term in the right hand side of (12) represent the loss of stems, either through natural loss or as a consequence of their sap having been sucked by larval A. itadori. The term is the rate of production of new shoots, and the last term represents the death, at time t, of stems which were new shoots at time and survived both natural loss and A. itadori activity over the time period . The exponential coefficients in that term represent the probabilities of surviving each of these forms of death over that time period.

Finally, we need an equation for the variable R(t) representing the knotweed rhizome biomass. One of the disadvantages of biocontrol using A. itadori is that the predator does not attack the root system which is strong and readily grows into new plants. Therefore, complete eradication cannot be expected. The root system is extensive can grow horizontally reaching up to 7 m from the parent plant, and be up to 3 m deep. Carbohydrates, made in the leaves, are carried back to the rhizome system which can store large quantities and are used for both top (stem) growth and root growth. We allow for loss of rhizome biomass at a per-capita rate , the subscript standing for rhizome or root. Given the resilience of the rhizome system of F. japonica, the value of is likely to be small. We let be the production rate of new rhizome biomass, which is taken to be a function of stem biomass S(t) since root growth depends on the production of carbohydrate by photosynthesis in the leaves, and leaf biomass is directly related to stem biomass. We propose the following simple equation for rhizome biomass R(t):

Model Analysis

The model to be solved consists of Eqs. (6), (7), (12) and (13) subject to the initial conditionsThe last two equations in (14) are constraints on initial data which reflect a compatibility with (5) and (11). Non-negativity for all time holds only for initial data satisfying the two constraints, but the use of non-negative initial data that does not satisfy the two constraints would have only a transient effect on the dynamics.

In view of the apparent capability of a knotweed stand to grow in an unlimited manner in the absence of any form of control, it seems not at all unreasonable to consider the situation when the birth rate functions for stems and rhizome are taken as linear, and so we take

Unbounded Knotweed Growth in Absence of A. itadori

If A. itadori is absent, then , the equations for the knotweed stem and rhizome biomass S(t) and R(t) reduce to a simpler form from which we can obtain the following single equation for R(t):where we have used the integral expression (11) for S(t), with . With the choices (15), Eq. (16) becomes, after a substitution,Solutions of the form exist whenever satisfies a characteristic equation. Since the delay term in (17) has a positive coefficient and only involves finite delay, we may apply the theory in Smith (1995) to conclude that the dominant eigenvalue is real. The characteristic equation isand we may consider only its real roots to determine whether R(t) grows or decays with time. The ratio in round brackets is decreasing in and simple calculus arguments show that the dominant real root of (18) is positive, so that the rhizome biomass R(t) grows exponentially, if . If the reversed inequality holds, then it decays exponentially. We have proved the following result.

Theorem 3.1

Suppose that there are no A. itadori and that and are chosen as in (15). Then, if the knotweed stand grows exponentially, while if then it decays exponentially.

From this very simple analysis, we may make some useful observations. The stand grows if is sufficiently large, as expected, and also if is sufficiently small (which is likely to be the case, since roots are hard to destroy). But the dependence on is more complicated. Note that, if is very small, then the condition for exponential growth is approximated by which only holds if is sufficiently large, i.e. stems live long enough before they die back. So, smallness of is not by itself sufficient for the plant’s survival (but smallness of is). It is hardly surprising that the plant has evolved to have such tough roots that are difficult to destroy.

Impossibility of Eradication When

We prove the following result which establishes that, if the knotweed stand grows exponentially in the absence of A. itadori, then it cannot be eradicated in the presence of A. itadori.

Theorem 3.2

Suppose that and are chosen as in (15), that is bounded and that . Then the knotweed infestation cannot be eradicated completely by A. itadori.

Proof

Suppose eradication occurs, then as . Immediately, we have from (13) with (15), from (5) using that is bounded, and from (7). However, by Theorem 3.1, in the absence of A. itadori, the knotweed biomass grows unboundedly under the condition , giving a contradiction.

Unbounded Knotweed Growth in Presence of A. itadori

In this subsection we again consider the case when the stem and rhizome production rates are linear and given by (15), but we now introduce biocontrol using the psyllid A. itadori whose birth rate function per stem need not be linear and in fact is kept as a general function satisfying certain properties. We will demonstrate that it is possible for a solution of (6), (7), (12) and (13) to have its four variables L(t), A(t), S(t) and R(t) growing exponentially with the same growth rate that is lower than in the case, considered in Sect. 3.1, when A. itadori is absent. The linear analysis we present here is of a slightly unusual character and is purely formal, but needs careful explanation. It is not linearised analysis about an equilibrium. Rather, it is a linear analysis in which we assume, a priori, that a solution exists in which all components have the same exponential growth rate. We exploit the fact that, for such a solution, as the variables L / S and A / S approach constants whose values are not immediately known but are determined later in terms of the growth rate . This approximation linearises the four equations even though the function is kept general. The result is a characteristic equation for having an unusual structure. The existence of a positive root under certain conditions is shown thereby confirming, a-posteriori, the existence of a solution with the supposed exponential growth.

It makes sense to assume the parameter values are such that the A. itadori population will grow in the situation when there are few adult A. itadori (i.e. A(t) is low) but unlimited stems (). In this situation, Eq. (7) approximates toThe solution A(t) will grow if the coefficient of the delay term exceeds that of the undelayed term (see Kuang 1993), thus we assume thatWe suppose the existence of a solution of (6), (7), (12) and (13) such thatwith and . Using expression (11) for S(t), and the linear stem and rhizome production rates (15), we can writeFor large times, and so (20) becomes, approximately,With , this becomesAlso, as the variables get very large, and and so the integral expression (5) for L(t) becomes, approximately,which, with and , becomesFinally, we need information about . This will be obtained from the equation for A(t), Eq. (7). As and so (7) becomes, approximately,With and , this becomeswhich implicitly determines . So, we let be defined as the solution (existence of which will be shown below) ofThen and is determined in terms of and using (24). Substituting the result into (22), we obtain the following characteristic equation to be solved for :with defined by (26). As regards existence of , we are only assured of this for small positive . Though zero is not an admissible value for in this subsection, existence of nevertheless follows from (19) and the other properties of . Continuity arguments then yield existence of for small positive . Moreover, is decreasing in for sufficiently small .

As a check on (27), we can let . Recall that h is the larva’s digestion time per unit biomass of sap consumed. Increasing this parameter should make the A. itadori larvae less effective as a biocontrol agent. Indeed, in the limit as , (27) formally tends to the earlier characteristic equation (18) for the case when A. itadori is absent.

Recall that the foregoing analysis presumes that is a small positive number. Continuity arguments yield the existence of a small positive root of (27) ifwhere means slightly less than, because if (28) holds then the graphs of the left- and right-hand sides of (27), as functions of , must cross each other at some small positive value of . Thus we have formally established the existence of a slowly growing exponential solution of (6), (7), (12) and (13) under condition (28). Since the analysis presumes that , it makes no sense to consider the existence of negative real roots of the characteristic equation (27), and no conclusions can be drawn from the existence thereof. However, in view of the impossibility of solutions tending to zero (Theorem 3.2), the nonexistence of positive roots of (27) could suggest that solutions remain bounded and this would be an important point in relation to the issue of effective biocontrol using A. itadori. We therefore conjecture that if exceeds the right-hand side of (28) then solutions remain bounded. In terms of the parameter h, which measures the insect’s digestion time per unit biomass consumed, intuitively we expect that effective biocontrol with bounded solutions will only happen for h below some threshold, with solutions growing exponentially otherwise. The function is decreasing in x, and therefore the right-hand side of (28) increases with h (note that does not involve h). This leads us to make the following conjecture.

Conjecture 3.3

If , where satisfiesthen solutions of (6), (7), (12) and (13) remain bounded.

The true threshold value of h distinguishing between boundedness and exponential growth turns out to be higher than , as we demonstrate later through numerical simulations. Again it is stressed that the foregoing analysis is purely formal.

Equilibrium Population for Small h

In this section, we assume parameters are such that the knotweed would grow exponentially in the absence of A. itadori. We show that, in the presence of A. itadori, if h is sufficiently small then solutions may evolve to a stable equilibrium so that the predator controls, but does not eradicate, the knotweed. No equilibrium exists if h is large. Therefore, A. itadori can only stabilize the knotweed population if it is able to digest the knotweed biomass sufficiently quickly. Otherwise, A. itadori can only slow down the rate of growth of the knotweed stand which still grows without bound. First note that, if (15) holds, the components of any equilibrium satisfyLet the function be defined as the positive root of the equationA suitable, and biologically realistic, assumption which guarantees that is well defined for each is (35) below:Assumption (35) basically states that the maturation rate per stem of A. itadori exceeds the death rate per stem if the ratio of adult insects to stems is low, but the opposite is true for high values of this ratio. The threshold in (35) is what we call . It is useful to keep in mind the realistic choice mentioned earlier, which satisfies (35) for large enough P.

Recall that the function decreases from at to at . If (35) holds, calculus arguments show that is well defined for any , and that increases with S and is confined between two finite positive values. Note also that . We prove the following result on the existence and stability of a positive equilibrium satisfying (30)–(33) for sufficiently small h.

Theorem 3.4

Suppose that (15) and (35) hold, and that (i.e. knotweed grows in the absence of A. itadori). Then, if h is sufficiently small, Eqs. (6), (7), (12) and (13) admit an equilibrium in which each component is strictly positive. Moreover, if and , where is defined by (34), then any real roots of the characteristic equation of the linearisation about this equilibrium must be negative. If h is too large then no such equilibrium exists.

Before we present the proof, we point out that the condition almost follows from (35) but is added to exclude the unlikely possibility of tangency. The condition states that the ratio must be on the increasing side of any non-monotone , while will hold if is concave () in the region of interest, which is realistic for the types of choices we have in mind.

Combining (30) and (31) and requiring givesLet be the root of the equationSince it is easily checked that exists and is unique. Note that and . Therefore, from (32) we obtain a single equation for the equilibrium component :The other equilibrium components are then determined in terms of fromThe existence or otherwise of a root of (38) can be considered from graphical arguments. Imagine the left- and right-hand sides are plotted against . If the left-hand side of (38), as a function of , decreases monotonically to a limit , while the right-hand side starts at 0 and grows without bound, becoming linear in for large. Therefore a root exists for . If is small, then the left-hand side initially decreases with but eventually starts to increase becoming approximately linear with a low growth rate of so that again (38) can be solved for a positive . However, if h is sufficiently large then the left-hand side increases for all , again becoming roughly linear but with a steeper growth rate than the right-hand side. This prevents the graphs from having an intersection and so there is no equilibrium.

To investigate the linear stability of the equilibrium that exists for small h, we set and . Collecting the linear parts of the resulting equations in terms of the new variables, and making use of the equilibrium equations (30)–(33), we obtainandSolutions are sought of the form . This leads, after some immensely complicated algebra, to a characteristic equation to be solved for . Lettingand casting the characteristic equation in a form where and have been eliminated, the result iswherewhere we have used (39) and we remind the reader that is the unique positive root of (37). Let us consider the behaviour of the above expressions for . Now is a positive decreasing function which satisfies . Also is decreasing so that . Moreover and . Using these facts, calculus arguments yield that the functions have the following properties for : and are positive and decreasing to zero as is increasing, as , and is increasing if h is sufficiently small.

We claim that the denominator of (45) is positive for all if h is sufficiently small. The first two terms in the right-hand side of (51) cancel out when , due to the fact that satisfies (37). Therefore, for small is a small negative number of O(h), while and are positive and O(1) in h. So for sufficiently small h. Also as , so the question is whether can go negative for some intermediate . Since is increasing, this evidently cannot happen after has become positive so we need only concern ourselves with a small compact interval near where . In every point of this interval, is small and of O(h), while and are positive and O(1) in h. Therefore, for sufficiently small throughout the small interval near where . It should be recalled that depends on h, however approaches a strictly positive value as .

We now have that, if h is sufficiently small, the right-hand side of (45) is negative for all . The left-hand side is increasing and is therefore positive for all , since by hypothesis. Therefore, any real roots of (45) must be negative and the proof is complete.

Numerical Simulations

In this section, we present the results of numerical simulations which confirm the analytical results of Sect. 3, and yield evidence for the truth of Conjecture 3.3.

We assume that the decreasing function is given by , that , and we use the following parameter values:and the following initial estimate of the number of adult A. itadori and quantity of rhizome biomass:For a complete list of the symbols and parameter values used, see Table 1. With the parameter values (52), we haveandwhere the unit of is the number of days per unit of stem biomass per larva.

If , we expect that the solutions remain bounded (Conjecture 3.3). Figure 1 shows that the total numbers of larval and adult A. itadori approach equilibrium values after a long oscillatory transient of around 3 years, with higher numbers of larval than adult A. itadori at equilibrium. The number of adults A(t) approaches a very small positive limit, not zero, and Fig. 1 therefore suggests that effective long-term control of F. japonica by A. itadori is possible with surprisingly few adult A. itadori. The oscillatory nature of the transient is not surprising since this is a predator–prey interaction between A. itadori larvae (the predator) and the F. japonica stems (the prey). The simulation shown in Fig. 2 is for a situation with but again indicates slow convergence to equilibrium, this time on an even longer timescale. Thus, is not a true threshold value for h; bounded or convergent behaviour is possible even for some (though not for all h). A detailed study of the transient region of Fig. 2 showed that the fluctuations are on a cycle of roughly 12 days, about the same as the mean adult longevity of adult A. itadori with the value for in (52). The bottom panel of Fig. 2 shows that both total stem biomass and total rhizome biomass approach equilibrium values after long transients. The variable S representing total stem biomass has a much higher amplitude of oscillation than the rhizome during the transient, no doubt because A. itadori eats only the stems. Mathematically, a lower amplitude of oscillation for rhizome biomass R can be anticipated because the solution R(t) of (13) at time t will involve a weighted average of past values of the stem biomass variable S.

Fig. 1

Solutions of model (6), (7), (12) and (13) subject to the initial conditions (14). The top panel shows the numbers of larval (left) and adult (right) A. itadori against time. The bottom panel shows total F. japonica stem (left) and rhizome (right) biomass against time. Parameter values used are given in (52) and . The variable A(t) tends to a small positive limit, not to zero, as

Fig. 2

Solutions of model (6), (7), (12) and (13) subject to the initial conditions (14). Parameter values used are given in (52) and . The variable A(t) tends to a small positive limit, not to zero, as

As we anticipated in the comment after Conjecture 3.3, which has been confirmed by the simulation shown in Fig. 2, when h is slightly above the total F. japonica stem and rhizome biomass variables remain bounded; in fact they continue to approach equilibrium values for h sufficiently close to . So the weed continues to be effectively controlled by A. itadori. Figure 3 shows a situation, for a larger , in which the variables still remain bounded but no longer approach equilibrium values. The situation here is one of sustained oscillations, the details of which are revealed in Fig. 4. If h is increased sufficiently, we eventually see unbounded oscillatory growth of F. japonica, as shown in Fig. 5, in which the peaks of all variables L(t), A(t), S(t) and R(t) grow in time. The total stem biomass still has a higher amplitude of oscillation than the rhizome. Such findings imply that the true threshold value for h, that distinguishes between boundedness and unbounded growth of solution variables, is higher than . This is not a surprise since most of the parameter values we used for F. japonica are based on approximations. All of these results together are very intuitive: if A. itadori voraciously consumes the stem biomass of F. japonica, then the growth of F. japonica is limited. Therefore, using A. itadori to control F. japonica is possible. However, F. japonica grows without bound if A. itadori is not sufficiently voracious.

Fig. 3

Solutions of model (6), (7), (12) and (13) subject to the initial conditions (14). Parameter values used are given in (52) and

Fig. 4

The situation of Fig. 3 showing the detail of the sustained oscillation

Fig. 5

Solutions of model (6), (7), (12) and (13) subject to the initial conditions (14). Parameter values used are given in (52) and

Discussion

We have developed a model to study the bio-control of F. japonica using one of its co-evolved natural enemies, the Japanese sap-sucking psyllid A. itadori. Our study focuses on a single isolated knotweed stand by modelling the growth of the stem and the rhizome with the stems under predation by larval A. itadori. Our final model consists of four differential equations with time delays subject to initial conditions with constraints. These four equations represent the rates of change of larval and adult A. itadori, and stem and rhizome biomass. Our results show that F. japonica will grow unboundedly in the absence of A. itadori if the natural loss of the stem biomass and the rhizome biomass is low. If such loss is high enough, then the knotweed stand decays exponentially. We have focused on the case when knotweed would grow in the absence of A. itadori, showing that in this case A. itadori cannot completely eradicate the knotweed but will either slow down its growth or, if it consumes the sap sufficiently voraciously, will cause the stem and rhizome biomass variables to remain bounded and possibly approach equilibrium values. The latter scenarios constitute the most desirable and effective form of biocontrol of knotweed whereby, based on a limited form of linearised analysis and numerical simulations, the knotweed and biocontrol agent coexist at a stable equilibrium or in the form of long-term sustained oscillations about an equilibrium.

It turns out that the dynamics of the model depend mainly on a parameter h, which measures how long it takes for an A. itadori larva to handle (digest) one unit of F. japonica stem biomass. Using model parameters based on the best available information, we provide an estimate of . If h is too large (i.e. h is much larger than ), then the model does not have a positive equilibrium and the plant biomass and insect numbers both grow together without bound, though at a lower rate than if A. itadori were absent. On the other hand, if h is sufficiently small, then the model implies that the knotweed plant biomass and insect numbers remain bounded and may approach equilibrium values. Rigorous mathematical analysis of our model is very challenging. In our study of the linearised stability of the coexistence equilibrium, we were only able to consider the real roots of the characteristic equation, even though the dominant root need not be real. In these circumstances, of course, negativity of the real roots does not conclusively yield stability of an equilibrium. However, the result on the existence of such an equilibrium for sufficiently small h is rigorous and this finding by itself is very important, since the knotweed stand grows without bound if no such equilibrium exists. It might be possible to prove boundedness of solution variables for sufficiently small h by other techniques, but our model is intractable to all strategies for proving boundedness that we have tried. Such a rigorous boundedness result for small h would be very desirable since the sole aim is to control the knotweed.

In this paper, we did not model the spatial spread of a knotweed stand explicitly, though we noted that the stand must spread in space if total stem biomass is increasing without bound with individual stems only growing to a certain height. Knotweed rhizomes are known to spread very vigorously and Smith et al. (2007) modelled the development of the rhizome system of a single stand using a 3D correlated random walk model. Their model predicts that the area of a stand increases quadratically with time and further suggests that it would be reasonable to model the growth of a stand using a reaction diffusion equation. Data from the spread of existing infestations could help to estimate the diffusivity. Development of our model to include the spread of a stand using such ideas is to be an area of future work. Spatial spread by other mechanisms such as downstream drift of root fragments (e.g., Hugh Dawson and Holland 1999) will also be important in future modelling efforts.