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 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 τ p (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 write which simply expresses the fact that larval A. itadori are lost at a per-capita rate μ l (S(t)) that depends on the total stem biomass. The function μ l (·) 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 μ ∞ l as its argument tends to infinity, that limit being the natural per-capita mortality when there is no shortage of sap. We also write Therefore l ξ (a) = l ξ (0) exp − a 0 μ l (S(η + ξ)) dη and, setting ξ = t − a, 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 b p (A(t)/S(t)) be the egg laying rate per stem which is a function of the number of adult A. itadori per stem. We will keep b p (·) as general as possible, but it is useful to keep in mind the frequently used prototype b p (x) = Pxe −qx where P, q > 0 (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 by If b p (·) is chosen as b p (x) = Pxe −qx , the implication is that the overall egg laying rate (i.e. for the whole knotweed stand) is l(t, 0) = P A(t)e −q A(t)/S(t) . This gives the expected formula P A(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 l(t, 0) = 0 if S(t) = 0 which again is correct because if there are no stems there are no leaves and no oviposition sites, and the factor e −q A(t)/S(t) measures the impact of crowdedness on egg laying. Using (3), (2) becomes and the total number of larvae is Differentiating this yields The 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 μ a 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 μ s . Letting s(t, a) be the age density of stem biomass, we may suppose that stem loss is governed by where 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 that Now, 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 b s (R(t)) be the production rate of stem biomass, where the subscript s stands for stems, and write Using (10) in (9), we find that the total stem biomass S(t) = τ s 0 s(t, a) da is, after a change in the variable of integration, given by where τ s 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 that The 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 b s (R(t)) 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 t − τ s and survived both natural loss and A. itadori activity over the time period [t − τ s , t]. 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 μ r , the subscript standing for rhizome or root. Given the resilience of the rhizome system of F. japonica, the value of μ r is likely to be small. We let b r (S(t)) 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 conditions The 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 A(t) = L(t) = 0, 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 L(t) = 0. With the choices (15), Eq. (16) becomes, after a substitution, Solutions of the form R(t) = exp(λt) 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 is and 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 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 b s (R) and b r (S) are chosen as in (15).
From this very simple analysis, we may make some useful observations. The stand grows if k r k s is sufficiently large, as expected, and also if μ r is sufficiently small (which is likely to be the case, since roots are hard to destroy). But the dependence on μ s is more complicated. Note that, if μ s is very small, then the condition for exponential growth is approximated by μ r < τ s k r k s which only holds if τ s is sufficiently large, i.e. stems live long enough before they die back. So, smallness of μ s is not by itself sufficient for the plant's survival (but smallness of μ r 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 b s (R) and b r (S) are chosen as in (15), that b p (·) is bounded and that μ r μ s < k r k s (1−e −μ s τ s ). Then the knotweed infestation cannot be eradicated completely by A. itadori.

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 b p (·) 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 t → ∞ 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 b p (·) 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 (S(t) → ∞). In this situation, Eq. (7) approximates to The solution A(t) will grow if the coefficient of the delay term exceeds that of the undelayed term (see Kuang 1993), thus we assume that We suppose the existence of a solution of (6), (7), (12) and (13) such that with λ > 0 and c i > 0, i = 1, . . . , 4. Using expression (11) for S(t), and the linear stem and rhizome production rates (15), we can write For large times, eσ L(t)/(1 + heσ S(t)) ∼ c 1 /(c 3 h) and so (20) becomes, approximately, With R(t) = c 4 e λt , this becomes Also, as the variables get very large, A(t)/S(t) ∼ c 2 /c 3 and μ l (S(t)) ∼ μ ∞ l and so the integral expression (5) for L(t) becomes, approximately, which, with L(t) = c 1 e λt and S(t) = c 3 e λt , becomes Finally, we need information about c 2 /c 3 . This will be obtained from the equation for A(t), Eq. (7). As t → ∞, A(t)/S(t) ∼ c 2 /c 3 and so (7) becomes, approximately, With A(t) = c 2 e λt and S(t) = c 3 e λt , this becomes which implicitly determines c 2 /c 3 . So, we let α(λ) > 0 be defined as the solution (existence of which will be shown below) of Then c 2 /c 3 = α(λ) and c 1 /c 3 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 α(λ) > 0, we are only assured of this for small positive λ. Though zero is not an admissible value for λ in this subsection, existence of α(0) nevertheless follows from (19) and the other properties of b p (·). Continuity arguments then yield existence of α(λ) > 0 for small positive λ. Moreover, α(λ) is decreasing in λ for sufficiently small λ. As a check on (27), we can let h → ∞. 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 h → ∞, (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) where < sl means slightly less than, because if (28) holds then the graphs of the leftand 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 λ > 0, 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 μ r 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 x → (1 − e −x )/x is decreasing in x, and therefore the right-hand side of (28) increases with h (note that α(0) does not involve h). This leads us to make the following conjecture.
The true threshold value of h distinguishing between boundedness and exponential growth turns out to be higher than h crit , 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 (L * , A * , S * , R * ) satisfy A suitable, and biologically realistic, assumption which guarantees that ϕ(S) is well defined for each S > 0 is (35) 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 ϕ(S). It is useful to keep in mind the realistic choice b p (x) = Pxe −qx mentioned earlier, which satisfies (35) for large enough P.
Recall that the function μ l (S) decreases from μ 0 l at S = 0 to μ ∞ l at S = ∞. If (35) holds, calculus arguments show that ϕ(S) is well defined for any S > 0, and that ϕ(S) increases with S and is confined between two finite positive values. Note also that ϕ(S * ) = A * /S * . We prove the following result on the existence and stability of a positive equilibrium satisfying (30)-(33) for sufficiently small h. (15) and (35) hold, and that μ r μ s < k r k s (1 − e −μ s τ s ) (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 (L * , A * , S * , R * ) in which each component is strictly positive. Moreover, if μ a > b p (ϕ(S * ))e −τ p μ l (S * ) , b p (ϕ(S * )) > 0 and b p (ϕ(S * )) − ϕ(S * )b p (ϕ(S * )) > 0, where ϕ(S) 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.

Theorem 3.4 Suppose that
Before we present the proof, we point out that the condition μ a > b p (ϕ(S * )) e −τ p μ l (S * ) almost follows from (35) but is added to exclude the unlikely possibility of tangency. The condition b p (ϕ(S * )) > 0 states that the ratio A * /S * must be on the increasing side of any non-monotone b p (·), while b p (ϕ(S * )) − ϕ(S * )b p (ϕ(S * )) > 0 will hold if b p (·) is concave (b p < 0) in the region of interest, which is realistic for the types of choices we have in mind.
Proof Combining (30) and (31) and requiring S * > 0 gives Let c = c * be the root of the equation Since μ r μ s < k r k s (1 − e −μ s τ s ) it is easily checked that c * > 0 exists and is unique. Note that eσ L * /(1 + heσ S * ) = c * and A * /S * = ϕ(S * ). Therefore, from (32) we obtain a single equation for the equilibrium component S * : The other equilibrium components are then determined in terms of S * from The existence or otherwise of a root S * > 0 of (38) can be considered from graphical arguments. Imagine the left-and right-hand sides are plotted against S * . If h = 0 the left-hand side of (38), as a function of S * > 0, decreases monotonically to a limit μ ∞ l c * /(eσ ), while the right-hand side starts at 0 and grows without bound, becoming linear in S * for S * large. Therefore a root S * > 0 exists for h = 0. If h > 0 is small, then the left-hand side initially decreases with S * but eventually starts to increase becoming approximately linear with a low growth rate of μ ∞ l hc * so that again (38) can be solved for a positive S * . However, if h is sufficiently large then the left-hand side increases for all S * > 0, 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 L = L * +L, A = A * +Ã, S = S * +S and R = R * +R. Collecting the linear parts of the resulting equations in terms of the new variables, and making use of the equilibrium equations (30)-(33), we obtaiñ Solutions are sought of the form (L,Ã,S,R) = (d 1 , d 2 , d 3 , d 4 )e λt . This leads, after some immensely complicated algebra, to a characteristic equation to be solved for λ. Letting and casting the characteristic equation in a form where L * , A * and R * have been eliminated, the result is where where we have used (39) and we remind the reader that c * is the unique positive root of (37). Let us consider the behaviour of the above expressions for λ ≥ 0. Now g(x) = (1 − e −x )/x is a positive decreasing function which satisfies g (x) > 0. Also μ l (·) is decreasing so that μ l (S * ) < 0. Moreover b p (ϕ(S * )) > 0 and b p (ϕ(S * )) − ϕ(S * )b p (ϕ(S * )) > 0. Using these facts, calculus arguments yield that the functions F 1 (λ), . . . , F 6 (λ) have the following properties for λ ≥ 0: F 1 (λ), F 2 (λ), F 4 (λ) and F 5 (λ) are positive and decreasing to zero as λ → ∞, F 3 (λ) is increasing, F 6 (λ) → 1 as λ → ∞, and is increasing if h is sufficiently small. We claim that the denominator F 2 (λ)F 5 (λ) + F 6 (λ) of (45) is positive for all λ ≥ 0 if h is sufficiently small. The first two terms in the right-hand side of (51) cancel out when λ = 0, due to the fact that c * satisfies (37). Therefore, for small h, F 6 (0) is a small negative number of O(h), while F 2 (0) and F 5 (0) are positive and O(1) in h. So F 2 (0)F 5 (0) + F 6 (0) > 0 for sufficiently small h. Also F 2 (λ)F 5 (λ) + F 6 (λ) → 1 as λ → ∞, so the question is whether F 2 (λ)F 5 (λ) + F 6 (λ) can go negative for some intermediate λ > 0. Since F 6 (λ) is increasing, this evidently cannot happen after F 6 (λ) has become positive so we need only concern ourselves with a small compact interval near λ = 0 where F 6 (λ) < 0. In every point of this interval, F 6 (λ) is small and of O(h), while F 2 (λ) and F 5 (λ) are positive and O(1) in h. Therefore, for sufficiently small h, F 2 (λ)F 5 (λ) + F 6 (λ) > 0 throughout the small interval near λ = 0 where F 6 (λ) < 0. It should be recalled that S * depends on h, however S * approaches a strictly positive value as h → 0.
We now have that, if h is sufficiently small, the right-hand side of (45) is negative for all λ ≥ 0. The left-hand side is increasing and is therefore positive for all λ ≥ 0, since F 3 (0) = μ a − b p (ϕ(S * ))e −τ p μ l (S * ) > 0 by hypothesis. Therefore, any real roots of (45) must be negative and the proof is complete.
We assume that the decreasing function μ l (S) is given by μ l (S) = k l e −S + m, that b p (x) = Pxe −qx , and we use the following parameter values: m = 0.0205, P = 16.9867, q = 1, μ a = 0.0806, e = 1, σ = 1, k s = 0.1, k r = 0.1, μ r = 0.01, μ s = 0.01, τ p = 32.2, τ s = 150, and the following initial estimate A 0 of the number of adult A. itadori and quantity R 0 of rhizome biomass: For a complete list of the symbols and parameter values used, see Table 1. With the parameter values (52), we have and where the unit of h crit is the number of days per unit of stem biomass per larva.
If h < h crit , 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 h > h crit but again indicates slow convergence to equilibrium, this time on an even longer timescale. Thus, h crit is not a true threshold value for h; bounded or convergent behaviour is possible even for some h > h crit (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 1/μ a of adult A. itadori with the value for μ a 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.
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 h crit the total F. japonica stem and rhizome biomass variables remain bounded; in fact they continue to approach  (52) and h = h crit − 1. The variable A(t) tends to a small positive limit, not to zero, as t → ∞ equilibrium values for h sufficiently close to h crit . So the weed continues to be effectively controlled by A. itadori. Figure 3 shows a situation, for a larger h > h crit , 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 h crit . 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  (14). Parameter values used are given in (52) and h = h crit + 5. The variable A(t) tends to a small positive limit, not to zero, as t → ∞ possible. However, F. japonica grows without bound if A. itadori is not sufficiently voracious.

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  (14). Parameter values used are given in (52) and h = h crit + 10 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 h crit . If h is too large (i.e. h is much larger than h crit ), 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.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
It is necessary to choose appropriate initial conditions which will be used as the history functions when one solves system (60) using dde23 in MATLAB. We choose the initial functions A 0 (θ ) and R 0 (θ ) to be constants A and R, respectively, and then use the integral constraints (67)-(70) to find constant initial data for the other variables, on their respective initial intervals, in terms of A and R. The constraints (67)-(70) determine that the constant initial functions for the variables L , S, X , and Y , which are used as the history functions in the MATLAB routine dde23, are the functions identically equal to the real numbers (again denoted L , S, X , and Y ) that satisfy X − exp − (k l e −S + m)τ p = 0, for specified values of A and R. Our simulations are of system (60), with initial data chosen to satisfy the integral constraints as just described.