A TOXIN-MEDIATED SIZE-STRUCTURED POPULATION MODEL : FINITE DIFFERENCE APPROXIMATION AND WELL-POSEDNESS

The question of the effects of environmental toxins on ecological communities is of great interest from both environmental and conservational points of view. Mathematical models have been applied increasingly to predict the effects of toxins on a variety of ecological processes. Motivated by the fact that individuals with different sizes may have different sensitivities to toxins, we develop a toxin-mediated size-structured model which is given by a system of first order fully nonlinear partial differential equations (PDEs). It is very possible that this work represents the first derivation of a PDE model in the area of ecotoxicology. To solve the model, an explicit finite difference approximation to this PDE system is developed. Existence-uniqueness of the weak solution to the model is established and convergence of the finite difference approximation to this unique solution is proved. Numerical examples are provided by numerically solving the PDE model using the finite difference scheme.

1. Introduction.The question of the effects of anthropogenic and natural environmental toxins on ecological communities is of great interest from both environmental and conservational points of view.Industrial toxins are one of the leading causes of pollution worldwide.Industrial toxins may arise as a result of air emissions, water releases, water seepage, air deposition or disposal and leaching of solid waste.The US Environmental Protection Agency has designated 126 priority pollutants [32] and the Canadian Council of Ministers of the Environment has a list of priority chemicals of concern for the protection of aquatic life [31].These priority substances can cause toxic effects when released into aquatic ecosystems.The effect of a toxic substance can, in principle, be exerted on all levels of the biological hierarchy, from cells to organs to organisms to populations to entire ecosystems.
To protect ecological environments and aquatic species, it is necessary to assess the risk of aquatic organisms exposed to toxins, and find relevant factors that determine the persistence and extirpation of organisms.Over the past several decades, 698 QIHUA HUANG AND HAO WANG mathematical models have been applied increasingly to predict the effects of toxins on a variety of ecological processes.These models include population models (single-population abundance, life history, individual-based, and metapopulation), ecosystem models (food-web, aquatic and terrestrial), landscape models, and toxicity-extrapolation models [10,12,24,25].The selection of specific models for addressing an ecological risk issue depends on the habitat, endpoints, and chemicals of interest, the balance between model complexity and availability of data, the degree of site specificity of available models, and the risk issue [24].A comprehensive review on the realism, relevance, and applicability of different types of models from the perspective of assessing risks posed by toxic chemicals is provided in [10,24].
Our search of the literature shows that toxin-dependent individual-based models and matrix population models are widely used to evaluate the ecological significance of observed or predicted effects of toxic chemicals on individual organisms and population dynamics.Despite the nonlinear dynamical nature of population-toxin interactions, relatively few differential equation models have been developed to describe population-toxin interactions (but see [11,13,14,15,19,20,23,28,29]).For those models that do exist, interactions are usually described by a system, which contains components representing the population density, the concentration of toxin in an organism, and the environmental concentration of toxin.Recently, we developed a toxin-dependent model given by a system of differential equations, to describe the impact of contaminants on fish population dynamics [16].Because the concentration of toxin in the environment is not affected significantly by mortality or metabolic processes of population, our toxin-dependent model focused on the impact of toxin on the population and ignores the influence of the population on the concentration of toxin in the environment.The concentration of toxin in the environment hence was treated as a parameter.The model was connected to literature-sourced experimental data via model parameterization of the toxic effects of methylmercury on rainbow trout (Oncorhynchus mykiss).The parameter estimates were then used to illustrate the long-term behavior of rainbow trout population.The numerical results provided threshold values of concentration of methylmercury in the environment to maintain populations and prevent extirpation.
It is significant that all above-mentioned population models assume that all individuals in a population have the same vital rates (reproduction, growth, and mortality rates) and the same sensitivity to toxins.However, in reality, different individuals may have different vital rates, and intake and egestion of toxins may depend on age, weight, and size.Motivated by these, in this work we extend the toxin-dependent ordinary differential equation model in [16] to a size-structured model with toxin effect.The model is given by a system of first order hyperbolic partial differential equations (PDEs) that includes two governing equations: one equation presents a generic description of the population growth under the influence of contaminant, while the other equation is the balance equation for the concentration of contaminant contained in individuals of the population.To our knowledge, this work presents the first derivation of a PDE model in the field of ecotoxicology.Although our toxin-mediated population model is developed in terms of size, the model and the results in this study are applicable to age-and weight-structured populations.
We are concerned with the existence and uniqueness of the weak solution of the system of PDEs.For this purpose, we develop an explicit finite difference approximation for the system in the spirit of the one initially used in [9,27] for conservation laws and later extended to nonlocal first order hyperbolic initial-boundary value problems arising in population ecology [1,5,6].In general, explicit schemes are computationally more practical and faster schemes for such problems (e.g., see [26]).Such a scheme results in an important tool: a numerical method that can be easily used for approximating solutions to this model, and for which convergence results are established.
The rest of the paper is organized as follows.In Section 2, we develop a toxinmediated aquatic population model.In Section 3, we define a weak solution of the system and develop an explicit finite difference approximation to the solution, we prove the convergence of approximation and existence-uniqueness of the weak solution.In section 4, we use the finite difference scheme to numerically solve the PDE model.Finally, a discussion section completes the paper.
2. Model formulation.Since we are interested an aquatic environment, the population is measured by concentration of biomass, rather than number.The state of the population is given by u(x, t), the size-distribution of the concentration of population biomass at time t.The formal meaning of u is that b a u(x, t)dx is the concentration of the biomass of individuals having size between a and b at time t.If there is no toxin, the population dynamics can be represented by the following McKendrick-vonFoerster equations [1,5,6,22]: with an appropriate initial condition.Here x min and x max denote the minimize size and maximum size of the population, respectively.The function g represents the growth rate of an individual of size x, µ denotes the mortality rate of an individual of size x, and β is the reproduction rate of an individual of size x.
Taking the effect of a toxin on the population into account, we let I(x, t) be the size-distribution of toxin in the population at time t.The formal meaning of I is that b a I(x, t)dx is the concentration of the toxin in the individuals having size between a and b at time t.Consider the rate of change of toxin in the population from size x to size x + ∆x, we have Rate of change of toxin in the individuals over the size interval (x, x + ∆x) = (rate of toxin entering the interval) − (rate of toxin leaving the interval) −(loss rate term)+(uptake rate term).This balance equation yields ( The growth of size leads to the individuals with certain concentration of the toxin enter and leave the small category.This leads to the terms I(x, t)g(x) and −I(x + ∆x, t)g(x + ∆x, t).The function E(t) is the concentration of toxin the environment at time t.The term aEu∆x represents the concentration of toxin uptaken by the population from the environment, which is modeled according to the Law of Mass Action, hence is proportional to both the concentration of toxin in the environment and the concentration of population biomass (see [16,28]), where a is the uptake coefficient.The term σI∆x represents the toxin elimination due to metabolic process of the population, where σ is per unit rate of toxin elimination.The term µI∆x represents the loss of the toxin due to death.Dividing (2) by ∆x and then letting ∆x → 0, we obtain Following [16,28], we consider the direct impact of the toxin on the population dynamics, which is realized by (body burden), given by I(x, t)/u(x, t) := v(x, t), concentration of toxin per unit population biomass.Using (1) and (3), we obtain the following equation for the body burden v: Let P (t) = xmax xmin u(x, t)dx be the total population biomass.We assume that the individual vital rates (reproduction, growth, and mortality rates) depend on the total population biomass due to competition for resources.Taking the effects of body burden on population vital rates into account, we consider the following nonlinear, nonlocal, hyperbolic PDE system that describes the dynamics of a population living in a polluted environment: u t + (g(x, P (t))u) x + µ(x, P (t), v(x, t))u = 0, (x, t) ∈ (x min , x max ) × (0, T ), In the absence of toxin, model (4) reduces to a classical size-structured population model of form (1), which has been well studied (see [1,6,22] and references therein).
3. Finite difference approximation and existence-uniqueness.In this section we establish the existence and uniqueness of weak solution to model (4).This is done through the following series of steps: 1) We construct a finite difference approximation for model (4).2) We establish some estimates for the solutions to the difference approximation.3) These estimates are then used to show that a set of functions generated from the difference approximation is compact in L 1 topology, which allows us to pass to the limit along a subsequence.This leads to the existence of a weak solution.4) Finally, we prove uniqueness of the weak solution, and hence establish convergence of the difference approximation.
3.1.Weak solution and finite difference approximation.Throughout the discussion we let and ω be a sufficiently large positive constant.We assume that the parameters in (4) satisfy the following assumptions: (A1) g : D 1 → R is a Lipschitz continuous function with Lipschitz constant L and satisfies sup (x,P )∈D1 g(x, P ) ≤ ω.Furthermore, g(x, P ) > 0 for x ∈ [x min , x max ) and g(x max , P ) = 0, and g x (x, P ) is Lipschitz continuous with respect to x and P for constant L.
Multiplying the first and second equations in (4) by ϕ(a, t) and ψ(x, t), respectively, and then formally integrating by parts and utilizing the initial and boundary conditions, we define a weak solution of (4) as follows: is called a weak solution to problem (4) if this set satisfies the following: We divide the intervals [x min , x max ] and [0, T ] into n and l subintervals, respectively.The following notation will be used throughout this paper: ∆x = (x max − x min )/n and ∆t = T /l denote the size and time mesh lengths, respectively.The mesh points are given by: We denote by u k j and v k j , the finite difference approximation of u(x j , t k ) and v(x j , t k ), respectively, and let We define the difference operators and the 1 and ∞ norms of u k and v k by We then discretize PDE system (4) using the following finite difference approximation with the initial conditions The following condition concerning ∆t and ∆x is imposed throughout the paper: (A9) Assume that ∆t and ∆x are chosen such that We can equivalently write (6) as the following system of linear equations: , from the first two equations of ( 7), one can easily see that under the assumption (A9), Thus, from the third equation of ( 7), we find u k+1 0 ≥ 0. That is, the difference system (7) has a unique solution satisfying From a biological point of view, it is very important that the numerical approximation preserves the nonnegativity of the solution.

3.2.
Estimates for the difference approximations.We first show that the difference approximation is bounded in 1 norm.Lemma 3.2.There exist positive constants M 1 and M 2 such that Proof.Multiplying the first equation of ( 7) by ∆x, summing over j = 1, 2, • • • , n, and noticing that g k n = 0, we have Using the boundary condition given in the third equation of (7) and assumption (A3), we get Treating the second equation of ( 7) similarly, and by the assumptions (A1) and (A5)-(A6), we find . which implies the second estimate.
We then establish ∞ bound on the difference approximation.
Lemma 3.3.The following estimates hold: , then from the third equation of ( 7) and the assumption (A3) we get Otherwise, suppose that for some 1 , then from the first equation of ( 7) and (A1) we have A combination of ( 8) and ( 9) then yields , then from the second equation of (7) and the assumptions (A5)-(A6), we get This leads to the desired result.
The next lemma shows that the approximations v k j has bounded total variation.Lemma 3.4.There exists a positive constant M 3 such that ) and apply the operator D − ∆x to the second equation of (7) to get By the assumption (A9), we have Multiplying the above equation by ∆x, and summing over the indices j = 2, 3, For j = 1, using the second and fourth equations of ( 7) and (A9) we have Adding ( 11) and ( 12) we get (13) Furthermore, we find By Lemma 3.3 and the assumptions (A4)-(A6), there exists a positive constant c 1 such that Applying the above inequality to (13), we arrive at which implies the estimate.
The next lemma is necessary to show that the approximations u k j have bounded total variation.Proof.We have from the second and third equations of (7) that Using g k n = 0 and the assumption (A3), we get Thus, by Lemmas 3.2-3.4and the assumptions (A1) and (A3), there exists a constant c 2 > 0 such that Moreover, using the assumption (A3), we find Hence, Furthermore, Thus, by Lemmas 3.2-3.3and the assumptions (A1)-(A2), there exists a constant Simple calculations yield Thus, by Lemmas 3.2 and 3.4 and the assumptions (A1) and (A4)-(A6), there exists a constant c 3 > 0 such that n j=1 Applying the bounds ( 16) and ( 17) to ( 15), we conclude that there exists a positive constant M 4 such that |u k+1 With the help of the above lemmas, we will show that approximations u k j have bounded total variation as well.The total variation bound plays an important role in establishing the sequential convergence of the difference approximation (6) to a weak solution of (5).Lemma 3.6.There exists a positive constant M 5 such that Proof.Set ξ k j = D − ∆x (u k j ) and apply the operator D − ∆x to the first equation of ( 7) to get By the assumption (A9), we have Multiplying the above equation by ∆x, and summing over the indices j = 2, 3, For j = 1, using the first equation of ( 7) and the assumption (A9), we have Adding ( 18) and ( 19), we get By Lemma 3.2, we get where xj ∈ [x j−1 , x j ], and xj−1 Therefore, by Lemmas 3.3-3.4,there exist positive constants c 4 and c 5 such that applying ( 21) and ( 22) to (20), we obtain The result now follows from the above inequality.
The next result shows that the difference approximations satisfy a Lipschitz-type condition in t.
Lemma 3.7.There exist positive constants M 6 and M 7 such that for any q > p, we have Proof.Summing the first equation in (7) over j and multiplying by ∆x, we obtain By Lemmas 3.2 and 3.6, and the assumptions (A1)-(A2), there exists a positive constant M 6 such that Thus, Similarly, by virtue of ( 17), we obtain

3.3.
Convergence of difference approximation and existence of a weak solution.Following [27] we define a family of functions {U ∆x,∆t } and {V ∆x,∆t } by Then by Lemmas 3.2-3.7 the set of functions ({U ∆x,∆t }, {V ∆x,∆t }) is compact in the topology of L 1 ((x min , x max ) × (0, T )) × L 1 ((x min , x max ) × (0, T )), respectively, and hence as in the proof of Lemma 16.7 on page 276 of [27], we have the following lemma.
Lemma 3.8.There exists a sequence of functions ) such that the limit functions satisfy The next theorem shows that the set of limit functions u(x, t), v(x, t) constructed via our difference scheme is actually a weak solution of problem (4).Theorem 3.9.The set of limit functions u(x, t) and v(x, t) defined in Lemma 3.8 is a weak solution of (4) and satisfies Proof.Let ϕ ∈ C 1 ((x min , x max ) × (0, T )) and denote the finite difference approximations ϕ(x j , t k ) by ϕ k j .Multiplying the first equation of the difference scheme (7) by ϕ k+1 j , we have Hence, Multiplying the above equation by ∆x, summing over k = 0, 1, • • • , l − 1 and j = 1, 2, • • • , n, and using g k n = 0 and the third equation of ( 7), we have (23) On the other hand, let ψ ∈ C 1 ((x min , x max )×(0, T )) and denote the finite difference approximations ψ(x j , t k ) by ψ k j .Multiply the second equation of ( 7) by ψ k+1 j to find Hence, Multiplying the above equation by ∆x , summing over k = 0, 1, • • • , l − 1 and j = 1, 2, • • • , n, and using g k n = 0 and v k 0 = 0, we have Using ( 23) and (24) and following an argument similar to that used in the proof of Lemma 16.9 on page 280 of [27] we obtain, by letting n, l → ∞, that the limit of the difference approximations defined in Lemma 3.8 is a weak solution of (4).Taking the limit in the bounds obtained in Lemmas 3.2-3.3,we get the bounds on P (t), u L ∞ ((xmin,xmax)×(0,T )) and v L ∞ ((xmin,xmax)×(0,T )) .
3.4.Uniqueness of the weak solution.The following theorem guarantees the continuous dependence of the solution u k j and v k j of ( 7) with respect to the initial condition u 0 j and v 0 j .
Theorem 3.10.Let {u k j , v k j } and {û k j , vk j } be the solutions of ( 7) corresponding to the initial conditions {u 0 j , v 0 j } and {û 0 j , v0 j }, respectively.Then there exists positive constants γ such that Then w k j , z k j satisfy the following: where ĝk i = g(x j , P k ) and similar notations are used for the rest of the parameters.Using the first equation of ( 25) and the assumption (A9), we obtain Multiplying the above inequality by ∆x, summing over the indices j = 2, 3, • • • , n and noticing that g k n = 0, we get For j = 1, by the second and fourth equations of ( 25) and the assumption (A9), we find Adding ( 26) and ( 27), we get Moreover, by the assumption (A3) , we have Note that Therefore, by the Lemma 3.3, there exists a positive constant c 6 such that Furthermore, where Therefore, by Lemmas 3.2-3.3and 3.6, and noticing that there exists a positive constant c 7 such that (31) Applying ( 29) and ( 31) to (28), we have On the other hand, using the second equation of ( 25) and the assumption (A9), we have Multiplying the above equation by ∆x, summing over the indices j = 1, 2, • • • , n, and noticing that g k n = 0 and z k 0 = 0, we have On the other hand, from (33)-(34), we find Adding ( 38) and (39), and letting , we obtain (37).Furthermore, (37) is equivalent to Hence, Now, from Theorem 3.9 we can take the limit in (40) to obtain where Â(•, t)} are the unique solutions of (36) with any set of given functions {P (t), B(t)} and { P (t), B(t)}, respectively.We then apply the estimate given in (41) for the corresponding solutions of (36) with two specific sets of functions {P (t), Q(t), B(t)} and { P (t), Q(t), B(t)} which are constructed using the limits obtained in Lemma 3.8 as follows: Thus, we get Figure 1.The three-dimensional dynamics of the solutions u(x, t) and v(x, t).Other parameters: a = 0.05, σ = 0.01, E = 2. biomass in the unstructured population model in [16].If data on the effects of a target toxin on the vital rates of a specific population is available, one can consider acute and chronic guideline developments by numerically solving model (4) using the finite difference scheme (7).We expect that our structured model developed in this work would produce more precise quantitative results than the unstructured model in [16]. 5. Discussion.Mathematical models are useful tools for evaluating the ecological significance of observed or predicted effects of toxic chemicals on individual organisms and population dynamics.Traditional ecotoxicological models assume that all individuals take up and excrete toxins in the same way and ignore the toxin transfer between individuals due to reproduction and growth.Motivated by the fact that depending on age, weight, and/or size, different individuals may have different sensitivities to toxins, we developed a toxin-mediated size-structured model that is given by a system of first order fully nonlinear hyperbolic PDEs.It is very likely that this model is the first derivation of PDE model in the area of ecotoxicology.Although our toxin-mediated model is based on an aquatic environment, the model and the results in this study are applicable to populations in terrestrial ecosystems.
To facilitate model analysis, we defined a weak solution and constructed an explicit finite difference approximation for the model.We proved the existenceuniqueness of the weak solution and the convergence of the finite difference approximation to the unique solution.This allows us to numerically solve the model using the finite difference scheme.The numerical solutions are then used to investigate the effects of toxins on the population dynamics.It is worth pointing out that the existing PDE solvers (such as pdepe and hpde in Matlab) can not be used to solve the nonlocal, fully nonlinear hyperbolic PDEs of form (4). Hence, the finite difference scheme we developed here provides an important numerical tool that can be used to approximate solutions of systems of nonlinear first order hyperbolic PDEs.
Model (4) assumed that the population vital rates g, µ, and β depend on P (t), the total population biomass at time t.If we assume that g, µ, and β are independent of P (t) and µ, β are independent of v, then model (4) reduces to a first-order linear system which can be solved using the method of characteristics.However, we are not interested in such a simplified version as it ignores the effects of toxin on population dynamics.If there is no toxin, model (4) reduces to a classical sizestructured population model which have been well studied in the literature.For example, a semilinear form of size-structured model (where g = g(x, t), β = β(x, t), and µ = µ(x, t, P (t))) was studied in [2,3].Therein, monotone approximations were developed, and upon passing to the limit a solution to the problem was obtained, whereas uniqueness was obtained via comparison results.For the quasilinear case (where g = g(x, P (t)), β = β(x, P (t)), and µ = µ(x, P (t))), the well-posedness was discussed in [8], wherein the method of characteristics together with a fixed point argument was employed to establish the existence and uniqueness of the solution to the model.It would be interesting to investigate whether or not above-mentioned methods are applicable to our toxin-mediated model (4).
Steps toward further research include the following: 1) Studying the asymptotic behavior of the model, that is, the existence and stability of equilibrium solutions, as shown by Figure 2. 2) Connecting the model to data by choosing a representative species and a representative toxin.In particular, applying the data on the effects of toxins on reproduction, growth, and survival of individuals to predict how different toxin levels affect the long-term behavior of the population.3) Our numerical examples assume that concentration of the toxin in the environment is a constant.In reality, the toxin concentration may vary over time due to a variety of factors, it would be interesting to study how time-dependent toxins affect the population dynamics.4) In practice, contaminant-induced changes in individuals' behavior may also lead to the change of population abundance.We expect that the main results we obtained in this study are robust, even though the details will certainly change if we include this factor in the model.5) Extending our single population model to interacting population model by considering the interactions between species such as competition, cooperation and predation [19,23].6) Stochastic models are generally more realistic than their deterministic counterparts, hence it might be worth extending model (4) to its stochastic counterparts to study the effects of environmental stochasticity [17,18,30] or demographic stochasticity [4] on population dynamics in a polluted environment.

Figure 2 .
Figure 2. Comparison of total population biomass, P (t) (top panel) and body burden of individuals having size 0.5, v(0.5, t) (bottom panel) between different toxin levels.Solid lines (E = 1), dashed lines (E = 4), dot lines (E = 10).Other parameters and initial conditions are the same as those in Figure 1.