Parameter estimation in a coupled system of nonlinear size-structured populations. Mathematical biosciences and engineering

A least squares technique is developed for identifying unknown parameters in a coupled system of nonlinear size-structured populations. Convergence results for the parameter estimation technique are established. Ample numerical simulations and statistical evidence are provided to demonstrate the feasibility of this approach.


Introduction
A typical direct problem for structured populations is to use the knowledge of underlying mechanism at individual level such as growth, mortality and reproduction rates to deduce the behavior at population level.This approach has been extensively studied for many kinds of models which include structured and non-structured populations.In practice, however, our knowledge of the vital rates may be incomplete [40].In fact, in many animal and plant populations the processes at the individual level are not accessible to direct observation [47].For example, for nonlinear structured models the dependence of reproduction and mortality rates on the total population is sometimes completely unknown [37].Even for linear structured models, one may not be able to obtain the exact dependence of the vital rates on the age or size structure [40].In these cases, one resorts to an inverse problem approach, namely to use knowledge about the behavior at the population level (e.g, observations of total population numbers) to deduce the underlying mechanisms at the individual level.
In recent years many researchers have focused their attention on developing methodologies for solving inverse problems governed by structured population models (e.g, [1]- [3], [12]- [17], [19]- [23], [25]- [34], [40]- [49]).In what follows, we briefly review some of the recent work on such inverse problems.For age-structured population models, several approaches have been developed to recover unknown individual vital rates.For example, in [40,43] a fixed point iterative technique was developed to determine the death rate from census data on the age distribution of the population.Therein, conditions on the data are given that lead to a unique solution.In [26] the authors formulated the inverse problem as an operator equation and the least squares method is then used to compute its solution.Due to the ill-posedness of the problem, a regularization technique was considered.In addition, the authors prove that the resulting scheme has a convergence rate of Hölder type.However, no numerical results were reported.A least squares approach was also adopted in [19] for a nonlinear agestructured population model to estimate unknown coefficients from a set of fully discrete observations of the population.Although the convergence of the computed minimizers to a minimizer of the least squares problem was established and numerical results were presented, for many real populations it is generally difficult to obtain discrete observations with respect to age, whereas other quantities such as total population number are easily obtained.In [25] a model describing the evolution in time of size/age structured population was considered.
A moving finite element method was used to study the identification problem for such a model.Convergence results for the parameter estimation technique were reported.In [30], by writing a linear age-structured model using the cumulative formulation approach (see e.g., [24]), the authors studied the inverse problem of identifying the birth and death rates from data on the total population size and the cumulative number of births.They also provided conditions on the data that guarantee the uniqueness of the solution to the inverse problem.
For size-structured population models, the least squares approach has been often used for parameter identification.For example, it was used in [15,16] to estimate the growth rate distribution in a linear size-structured population model.A similar technique was subsequently applied to a semilinear size-structured model in [34] where the mortality rate depends on the total population due to competition.In [2] an inverse problem governed by a phytoplankton aggregation model was studied.Convergence and numerical results for identifying the coagulation kernel were provided.Later, this technique was extended to identify parameters in a size-structured population model in [1,3] where all the individual vital rates (growth, mortality and reproduction) depend on the total population level.Therein, these parameters are identified from a set of observations corresponding to the total population number.A finite difference method was then used to approximate the infinite dimensional problem.Convergence results for the computed parameter estimates to the true parameter were established.To our knowledge, [3] was the first paper to provide convergence results for parameter estimates when the growth rate is a nonlinear function of the total population (i.e., the size-structured model is represented by a quasilinear first order hyperbolic initial boundary value problem).
In this paper we extend the discussion in [3] to the following coupled system of quasilinear size-structured populations model: t + (g I (x, P (t; q))u I ) x + m I (x, P (t; q))u I = 0, (x, t) ∈ (0, L] × (0, T ], g I (0, P (t; q))u I (0, t; q) = C I (t) + N J=1 L 0 γ I,J β J (x, P (t; q))u J (x, t; q)dx, t ∈ (0, T ], (1.1) Here q = (q 1 , q 2 , . . ., q N ) with q I = (g I , m I , β I , C I ), I = 1, 2, . . ., N , the parameters to be identified.The function u I (x, t; q), I = 1, 2, . . ., N , is the parameter-dependent size density (number per unit size) of individuals in the Ith population having size x at time t, and The model (1.1), which was developed by the authors in [4], is a generalization of several size-structured population models (usually referred to as structured models with rate distributions) which have been investigated in [14,15,16,34].Motivated by the fact that, in addition to observable characteristics such as age or size of the individuals, non-observable genetic characteristics may often play a crucial role in the development of the individuals, researchers in [14] presented the first such generalization of the classical Sinko-Streifer model.
This model, which is a linear version of (1.1), has vital individual rates that are independent of the total population and distributed over an an infinite-dimensional admissible parameter space with a probability measure.It was shown through numerical simulations in [14] that there is a crucial difference between the dynamics of distributed rate size-structured population models and the classical Sinko-Streifer models.In particular, the classical Sinko-Streifer model cannot have dispersion of the density of the population in age or size except under biologically unreasonable conditions on the growth rate [15].That is why the classical Sinko-Streifer models are in conflict with field data collected by experimental biologists.These data sets show that a population with unimodal distribution evolves into a bimodal distribution (see [14] and [41]).In [17] the authors used least squares approach to fit these distributed rate models to data obtained in [14].The resulting good fit indicates that the need for such modification is crucial if these models were to be used as prediction tools.
In addition to extending the theory in [3] to the coupled quasilinear system (1.1), a main novelty of our current research is that we report on extensive numerical simulations.These simulations are then used to obtain statistical results (in the form of confidence intervals) which provide solid evidence on the feasibility of this approach.It is worth pointing out that with the exception of [28] the above-mentioned articles do not report on any statistical studies.
As the use of numerical methods for estimating functional parameters becomes more widely accepted in the biological sciences, it is becoming increasingly important for investigators to support the efficacy of proposed numerical algorithms with not only numerical simulation results but also confidence intervals on estimated parameters.This can be done by calculating standard errors in a number of sophisticated ways (e.g., pointwise confidence intervals or bands as in [38,39,48], uniform bands [32], simultaneous confidence bands [31], etc.).Here we simply compute the pointwise standard errors using the pointwise sample variances from a large (1000) number of inverse problem simulations.While in our efforts we emphasize (regularized) ordinary least square estimators, the ideas and methods presented in this paper can readily be used with maximum likelihood estimators as well as other standard estimators found in the statistical literature.
It is also worth noting another connection between statistical methods and our efforts in this paper.The models we use here involve a form of "mixing" distributions found in the literature on mixed effects, random effects or hierarchial methods (see for example, [20,21,22,35,36,46]). However the models we investigate entail mixing that cannot be decoupled into individual dynamics and thus result in fully coupled dynamics (see our closing remarks in Section 4).
By a weak solution to problem (1.1) we mean a bounded and measurable function u(x, t; q) = (u 1 (x, t; q), u 2 (x, t; q), . . ., u N (x, t; q)) satisfying γ I,J β J (x, P (s; q))u J (x, s; q)dx ds for t ∈ [0, T ], I = 1, 2, . . ., N , and every test function We first impose a condition on the initial data: for any I = 1, 2, . . ., N (H1) u I,0 ∈ BV [0, L] and u I,0 (x) ≥ 0. where Depending on the values of the constants 0 ≤ γ I,J ≤ 1, the model (1.1) may have two different interpretations.If γ I,I = 1 and γ I,J = 0, I = J, the model represents the dynamics of several populations competing for common resources.On the other hand, if γ I,J > 0, I, J = 1, 2, . . ., N , then the model may describe the dynamics of one population consisting of N subpopulations, each with its own characteristics.Hence, γ I,J represents the probability that an individual of the Jth subpopulation will reproduce an individual of the Ith subpopulation.Therefore, two different ways for observing data will be considered.
These lead to the following two different least-squares functionals to be minimized: The first one is based on the assumption that the model (1.1) describes N different competing populations.Hence observations Z I,k which correspond to the total number of individuals in the Ith population at time t k are assumed to be available (this case corresponds to γ I,I = 1 and γ I,J = 0, I = J).We define the least-squares cost functional for this case to be which is minimized over Q.The other case assumes that (1.1) models one species which has been divided into N not readily distinguishable subpopulations.In this case, we assume that we can only observe aggregate data Z k , the total number of individuals at time t k (this case corresponds to γ I,J > 0, I, J = 1, 2, . . ., N ).We define the least-squares cost functional which is minimized over Q.We remark that minimizing (1.4) over Q is equivalent to the maximum likelihood esti- are i.i.d.normal, and minimizing (1.5) over Q is equivalent to the maximum likelihood estimation of q if k = log The paper is organized as follows.In Section 2, we present a finite difference scheme for computing the solution of (1.1) and then provide convergence results for the parameter estimation technique.In Section 3, we give ample numerical and statistical results.Some concluding remarks are made in Section 4.

Approximation Scheme and Convergence Theory
The following notation will be used throughout the paper: ∆x = L/n and ∆t = T /l denote the spatial and time mesh size, respectively.The mesh points are given by x j = j∆x, j = 0, 1, 2, . . ., n and t k = k∆t, k = 0, 1, 2, . . ., l.We denote by u I,k j (q) and P k (q) the finite difference approximation of u I (x j , t k ; q) and P (t k ; q), respectively, and we let We define the difference operator and the 1 , ∞ and the BV norms of u I,k by We then discretize the partial differential equation in (1.1) using the following implicit finite difference approximation with the initial condition If we define then (2.1) can be equivalently written as the following system of linear equations for where T and A k is the following block diagonal matrix: with the lower triangular matrix Note that using the assumptions on our parameters one can easily show that equation (2.2) has a unique solution satisfying u k+1 (q) ≥ 0, k = 0, 1, . . ., l − 1.
The above approximation can be extended to a family of functions {U I ∆x,∆t (x, t; q)} defined by Since our parameter set is infinite dimensional, a finite dimensional approximation of the parameter space is also necessary for computing minimizers.To this end, we consider the following finite-dimensional approximations of (1.4) and (1.5), respectively: and each of which is minimized over Q M , a compact finite-dimensional approximation of the parameter space Q.In order to establish the convergence results for the parameter estimation technique, we use a similar approach to that in [3], which is based on the abstract theory in [18].
Theorem 2.1 Let q i = (q 1,i , q 2,i , . . ., q N,i ) and suppose that for each I, q I,i → q I in Q I and denote the solution of the finite difference scheme, and let u(x, t; q) = (u 1 (x, t; q), u 2 (x, t; q), . . ., u N (x, t; q)) be the unique weak solution of our problem with initial condition and parameter q, then U From the fact that Q I is compact and the results of [4], there exist positive constants c 1 , c 2 , c 3 , c 4 such that for each I = 1, 2, . . ., N , we have , where r > s.Thus, for each I there exists a BV ([0, L] × [0, T ]) function ûI (x, t) such that U I ∆x i ,∆t i (x, t; q i ) → ûI (x, t) in L 1 (0, L) uniformly in t.Hence, from the uniqueness of bounded variation weak solutions stated in [4], we only need to show that û(x, t) = (û 1 (x, t), û2 (x, t), . . ., ûN (x, t)) is the weak solution corresponding to the parameter q.To this end, we multiply the first equation of (2.1) by Multiplying the above equality both sides by ∆x i ∆t i and summing over j = 1, 2, . . ., n, k = 0, 1, . . ., l − 1, we find Since g I,k,i n = 0 and q I,i → q I as i → ∞ in Q I , passing to the limit we have Thus, û(x, t) is the weak solution corresponding to the parameter q.
Since the logarithm function is continuous on [1, ∞), as an immediate consequence of Theorem 2.1, we obtain the following: Corollary 2.2 Let U ∆x,∆t denote the numerical solution of (2.1) with parameter q i → q and ∆x i , ∆t i → 0. Then In the next theorem, we establish the continuity of the approximate cost functional, so that the computational problem of finding approximate minimizer is well-posed.
Theorem 2.3 Let ∆x and ∆t be fixed.For each q I ∈ Q I , let U I ∆x,∆t (x, t; q) denote the solution of the finite difference scheme, and q I,i → q Proof.Define {u I,k,i j } and {u I,k j } to be the solution of the finite difference scheme with parameter q i and q, respectively.Let v I,k,i j = u I,k,i j − u I,k j , then v I,k,i j satisfies the following: for 1 ≤ j ≤ n, and where P k,i denotes P k (q i ).Multiplying both sides of (2.6) by sgn(v I,k+1,i j )∆x and summing over j = 1, 2, . . ., n, we obtain (2.8) Using the fact for any a j with a j ≥ 0, j = 0, 1, 2, . . ., n, we have (2.9) By (2.7), we have (2.10) Summing (2.8) over I = 1, 2, . . ., N , and using (2.9) and (2.10), we obtain Noticing that g I,i (x j , P k,i ) − g I (x j , P k ) ≤ g I,i (x j , P k,i ) − g I,i (x j , P k ) + g I,i (x j , P k ) − g I (x j , P k ) , we have from (H4) the following: Similarly, we can show that max Furthermore, straightforward computations yield Hence, from (H4) we obtain max Then, we have Since for each k, ρ k,i → 0 as i → ∞, the desired result easily follows from this inequality.
Theorem 2.4 Suppose that Q M is a sequence of compact subsets of Q.Moreover, assume that for each q ∈ Q, there exists a sequence of q M ∈ Q M such that q M → q as M → ∞.Then the functional J ∆x,∆t has a minimizer over Q M .Furthermore, if q i M denotes a minimizer of J ∆x i ,∆t i over Q M and ∆x i , ∆t i → 0, then any subsequence of q i M has a further subsequence which converges to a minimizer of J .
Proof.The proof of this theorem is a direct application of the abstract theory in [18], based on the convergence of J ∆x i ,∆t i (q i ) → J (q).

Numerical Results
In this section, we present ample numerical simulations and statistical results.In all of the simulations below we assume L = 1, T = 1, and C I (t) = 0 for I = 1, 2, . . ., N .
In subsections 3.1 and 3.2, we assume N = 1 and that all the parameters are known except for β.To estimate β we use data which are generated computationally as follows: Let and we solve (2.1) and (2.3) for U ∆x,∆t (x, t).We set the data where ε k is a random sample from a normal random number generator with mean zero and standard deviation σ = 0.02.
3.1 1 − D linear estimation problem for finite dimensional parameter space when N = 1 In our first example we assume that β is of a separable form given by β(x, P ) = b(x) exp(−3P ), where b(x) = µx(1 − x ν ) with µ and ν two unknown constants to be identified.Hence, the solution to our least-squares problems involves identifying the two constants µ and ν from a compact subset of R 2 + so as to minimize the least-squares cost functional In order to test the performance of the parameter-estimation technique when no infinite dimensional effects are present, in Figure 1 we choose ∆x = ∆t = 0.005 for both generating the data and the numerical solution (2.3) in the least-squares problem.This avoids the infinite-dimensional effect of the partial differential equation given in (1.1).In fact, if the noise is removed from the data, and the parameters µ and ν are known, then numerically solving our model produces the exact data.
In Figure 2 we use ∆x = ∆t = 0.005 to generate the data while we use ∆x = ∆t = 0.01 for the numerical solution (2.3) in the least-squares problem.Thus, in this case the data are not exactly attained by our model even if the noise is removed (an error is present due to the finite-dimensional approximation of our infinite-dimensional model).The results of Figure 2 are obtained by using the same values for the rest of the parameters as those of Figure 1.
A similar format for presenting the results of 1000 inverse problem calculations was used in Figure 1 and 2. The left part of each of the figures represents the S (for our case S = 1000) numerical results for the estimated parameter b s (x) (s = 1, 2, . . ., S) versus the exact b(x), where these 1000 distinct numerical results graphed were obtained by solving 1000 inverse problems, each of which corresponds to a given noise sample { k }.The right part represents the figure of the corresponding 95% confidence interval (dashed line) versus the exact b(x) (solid line), where the 95% confidence interval is obtained by choosing the band between the upper 2.5% and lower 2.5% of these 1000 numerical results.Table 1 provides denotes the sampling standard error at the point x.Note that this is simply the usual asymptotic formula for the pointwise standard error (e.g., see p. 28, 37 of [21] and p. 308 of [45]).
Although the estimates in both figures are good, the results in Figures 1-2 and Table 1 suggest that infinite-dimensional effects can lead to a slightly under biased estimator.
We suspect that this bias depends on the choice of the numerical scheme used for solving the infinite-dimensional partial differential equation model.Here we are using an upwind scheme for approximating the model and a right-hand sum for approximating all the integrals involved.This biased estimator may be improved if, for example, a centered finite difference approximation is used together with a trapezoidal rule for integration.The above statistical results (essentially on how measurement error affects estimates) are based on a large number of numerical simulations (somewhat in the spirit of Bayesian based MCMC calculations used to estimate means and variances in a probability distribution from "experimental" data).Any estimate of model parameters from data can also be accompanied by an estimate of uncertainty using standard regression formulations from statistics [21].Thus, in the remaining part of this subsection, we present a statistical based method to actually compute the variance in the estimated model parameters q = (µ, ν).x AB(x) RAB(x) SE(x) 0.1 -0.0037 -0.6870 0.0749 0.2 -0.0092 -0.9580 0.0993 0.3 -0.0107 -0.8463 0.0975 0.4 -0.0079 -0.5497 0.0860 0.  1 and Figure 2, respectively.
Substituting P (t i ; q), P µ (t i ; q) and P ν (t i , q) for P (t i ; q), P µ (t i ; q) and P ν (t i ; q) in (3.1), respectively, we obtain the following approximation of X(q): 1 + P (t 1 ; q) P ν (t 1 ; q) 1 + P (t 1 ; q) P µ (t 2 ; q) 1 + P (t 2 ; q) P ν (t 2 ; q) Under standard assumptions of classical nonlinear regression theory, we know that if ˆ i ∼ N (0, σ 2 ), where ˆ i is the difference between observation and model at time t i , then the least-squares estimate q * is expected to be asymptotically normally distributed.In particular, for large samples, we may assume where q 0 is the true vector of parameters and σ 2 {X T (q 0 )X(q 0 )} −1 is the true covariance matrix (see [21], Chapter 2).
Since q 0 and σ 2 are not available, we follow a standard statistical practice [5]: substitute the computed estimate q * for q 0 and approximate 2) to obtain the standard deviation for our estimates.In particular, if then we take V 11 and V 22 to be the standard deviation for parameters µ and ν, respectively.The following two tables are the standard deviation of µ and ν for the results of the first eight numerical simulations of Figure 1 and Figure 2, respectively.2.
Table 4 provides the average standard deviation of µ and ν for the results of all the 1000 numerical simulations of Figure 1 and Figure 2, respectively.We note that in most practical situations using experimental data, one does not expect to have 1000 experiments performed.But the above procedures will produce estimates of variances even in the case when one has only one data set!3.2 1 − D linear estimation problem for infinite dimensional parameter space when N = 1 In this example, we assume that β is of a separable form given by β(x, P ) = b(x) exp(−3P ), where b(x) is an unknown parameter that we want to identify. Let Choose the parameter space Q = D. Clearly, by Arzela-Ascoli Theorem [33] Q is compact in C[0, 1].We approximate the infinite dimensional parameter space as follows: For M a positive integer and b ∈ Q, we set where φ i M (x; 0, 1) are the linear spline functions on a uniform mesh of the interval [0, 1].These are defined by . It can be readily argued that lim then the solution of our finite dimensional identification problem involves identifying the + so as to minimize the leastsquares cost functional (2.4).
In order to indirectly implement the compactness constraints of Q, we use a regularized least squares cost functional of the form where α > 0 is the regularization parameter.
The left part of each of the following figures again represents the S (=1000) numerical results of the estimated parameter versus the exact parameter b(x).The right part represents the figure of the corresponding 95% confidence interval (dashed line) versus the exact b(x) (solid line).The tables provide statistical results for the corresponding graphs.

Effect of infinite-dimensional model on parameter estimate.
In Figure 3 we use ∆x = 0.005 and ∆t = 0.005 to generate the data and the numerical solution (2.3) for the least-squares problem.This removes the infinite-dimensional effect of the partial differential equation given by (1.1).However, in Figure 4 we use ∆x = ∆t = 0.005 to generate the data and ∆x = ∆t = 0.01 to compute (2.3).Thus, in this case the data are not exactly attained by our model even if the noise is removed.We observe that while the estimates in both figures are good, the results in Figures 3-4 and Table 5 suggest that infinite-dimensional effects can lead to a slightly under biased estimator.x AB(x) RAB(x) SE(x) 0.  3 and Figure 4, respectively.
Effect of regularization parameter α on parameter estimate.In Figures 5 and 6 we change the parameter α while keeping the rest fixed.Clearly, low regularization parameter leads to relatively bad estimates although the estimator in this case seems to be the least biased (see Figure 5 and left part of Table 6).Increasing the value of α leads to better parameter estimates, but the estimator becomes more under biased (see Figure 6 and right part of Table 6).If this value is increased more, the estimator is more biased.Also the parameter estimate becomes worse than before.This suggests, not surprisingly, that there is an optimal choice for the parameter α which produces the best results for the parameter estimates.5 and Figure 6, respectively.In this section, we assume N = 2 and that all the parameters are known except for β 1 and β 2 .To estimate β 1 and β 2 , we assume that they are of a separable form given by β 1 (x, P ) = b 1 (x) exp(−P ) and β 2 (x, P ) = b 2 (x) exp(−P ), respectively, where b 1 (x) and b 2 (x) are unknown parameters to be identified.To estimate b 1 (x) and b 2 (x), we use data which are generated computationally as follows: Let γ I,J = 1, I = J 0, I = J for Figure 7 and γ I,J = 0.5, I, J = 1, 2 for Figure 8, u I,0 (x) = 3 exp(−2(x − 0.1) 2 ), and for the parameters g I , m I and β I we use the following choice of functions: and solve (1.1) for U I ∆x,∆t (x, t), I = 1, 2. We set the data 7 and 8, where ε I,k and ε k both are the random sample from a normal random number generator with mean zero and standard deviation σ = 0.02.
We choose the parameter space We approximate the infinite dimensional parameter space as follows: For M 1 , M 2 positive integers and any (b Clearly, lim then the solution of our finite dimensional identification problem involves identifying the + so as to minimize the least-squares cost functional (2.4) or (2.5).
In order to indirectly implement the compactness constraints of Q, we still use the regu-larized least-squares cost functional.For Figure 7 we use the form and for Figure 8 we use the form where α I > 0, I = 1, 2 are the regularization parameters and m = 100 for Figures 7 and 8.
In the rest of our simulations we use ∆x = ∆t = 0.005 to generate the data and ∆x = ∆t = 0.01 to solve the least-squares.Thus, in these cases the data are not exactly attained by our model even if the noise is removed.Note that the results in Figure 7 and Table 7 are slightly better than those in Figure 8 and Table 8.This is expected since in Figure 7 we are sampling data for each of the two populations, which provides more information than sampling the sum of the two populations only, as is the case in Figure 8. Also note that in both of these figures we let M = M 1 = M 2 = 10.

Concluding Remarks
In this paper we have developed a numerical technique for identifying unknown parameters in a general size-structured population model.A main focus of the paper is on a statistical study of the parameter estimation technique.This was carried out by calculating pointwise standard errors on the estimated parameters (functions) via use of thousands of numerical experiments.
Several conclusions can be drawn from our studies.1) The method discussed above seems to perform well and produce good confidence intervals for the parameters.2) When the infinite dimensional effects of the model and the parameter space are removed, the resulting numerical and statistical values suggest that the least-squares technique produces very good unbiased parameter estimates.3) The type of numerical scheme used for approximating the infinite-dimensional model as well as the parameter space may influence the bias in the parameter estimation technique.4) The commonly used regularization term is crucial for enforcing compactness and obtaining better estimates.However, it may also introduce more bias in the estimator.
We note in closing that the system (1.1) investigated in this paper is a special case of the measure dependent aggregate dynamics problems formulated in [6] wherein individual (uncoupled) dynamics are not available.Inverse problems for such systems have been investigated in a number of applications including cellular level HIV modelling [7], hysteresis in viscoelastic materials [8,9], shear waves in biotissue [10], and electromagnetic interrogation in complex materials [11].In a more general formulation (currently under investigation by the authors), one has a probability distribution F of individual parameters q(x, P ) = q = (g, m, β, C) on an admissible set Q.The system (1.1) is replaced by a continuum of systems for u(x, t; q(x, P )) with the total population P (t; F ) given by P (t; F ) = Q L 0 u(x, t; q)dx dF (q) = Q L 0 u(x, t; q)dx f (q)dq, the latter equality holding if F has a density f .The aggregate dynamics for u depend explicitly on F through the dependence of the individual rate parameters (g, m, β, C) on the total population P .
If F is a discrete measure with N atoms at q J of mass f J , then we have f J L 0 u(x, t; q J )dx.

2 )
is the total population at time t.The function g I denotes the growth rate of an individual in the Ith population, m I denotes the mortality rate of an individual in the Ith population, and β I is the reproduction rate of an individual in the Ith population.The function C I represents the inflow rate of the Ith population of zero-size individuals from an external source (e.g., in a tree population model seeds moved by wind).
statistical results for the corresponding graphs, where AB(x) = 1 S S s=1 (b s (x)−b(x)) denotes the average bias for all approximations at x, RAB(x) = 100 AB(x) b(x) denotes the relative average bias for all approximations at x and

Figure 1 :
Figure 1: ∆x = ∆t = 0.005 to generate the data and solve the least-squares.For the left part of the figure, each of the grey lines (....) denotes a distinct result for a given sample { k }.

Figure 2 :
Figure 2: ∆x = ∆t = 0.005 to generate the data and ∆x = ∆t = 0.01 to solve the leastsquares.For the left part of the figure, each of the grey lines (....) denotes a distinct result for a given sample { k }.

Figure 1 Table 4 :
Figure 1 Figure 2 µ 1.1921 1.9197 ν 0.4566 0.8572Table 4: Average of standard deviation for all the results of the numerical simulations of Figures 1-2.

Figure 3 :
Figure 3: M = 10, α = 3e − 5.Each of the grey lines (....) of the left part of the figure denotes a distinct result for a given sample { k }.

Figure 4 :
Figure 4: M = 10, α = 3e − 5.Each of the grey lines (....) of the left part of the figure denotes a distinct result for a given sample { k }.

Figure 5 :
Figure 5: M = 10, α = 1e − 5.Each of the grey lines (....) of the left part of the figure denotes a distinct result for a given sample { k }.

Figure 6 :
Figure 6: M = 10, α = 5e − 5.Each of the grey lines (....) of the left part of the figure denotes a distinct result for a given sample { k }.

3.3 1 −
D linear estimation problem for infinite dimensional parameter space when N = 2 The upper-left part and the lower-left part of the following two figures represent the S (=1000) numerical results of the estimated parameters b 1 M 1 (x) and b 2 M 2 (x) versus the exact parameters b 1 (x) and b 2 (x), respectively.The upper-right part and the lower right part represent the figures of the corresponding 95% confidence interval (dashed line) versus the exact b 1 (x) and b 2 (x) (solid line), respectively.The tables provide statistical results for the corresponding graphs.

Table 1 :
Left and right tables are statistical results for Figure

Table 2 :
Standard deviation for the results of the first 8 numerical simulations of Figure1.

Table 3 :
Standard deviation for the results of the first 8 numerical simulations of Figure

Table 5 :
Left and right tables are statistical results for Figure

Table 6 :
Left and right tables are statistical results for Figure