Sixth Mississippi State Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, Conference 15 (2007), pp. 67–75. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu ftp ejde.math.txstate.edu (login: ftp)

NONSTANDARD NUMERICAL METHODS FOR A CLASS OF

PREDATOR-PREY MODELS WITH PREDATOR INTERFERENCE

DOBROMIR T. DIMITROV, HRISTO V. KOJOUHAROV

Abstract. We analyze a class of predator-prey models with Beddington-DeAngelis type functional response. The models incorporate the mutual interference between predators, which stabilizes predator-prey interactions even when only a linear intrinsic growth rate of the prey population is considered. Positive and elementary stable nonstandard (PESN) finite-difference methods, having the same qualitative features as the corresponding continuous predatorprey models, are formulated and analyzed. The proposed numerical techniques are based on a nonlocal modelling of the growth-rate function and a nonstandard discretization of the time derivative. This approach leads to significant qualitative improvements in the behavior of the numerical solution. Applications of the PESN methods to specific Beddington-DeAngelis predator-prey systems are also presented.

1. Introduction

The traditional mathematical model of predator-prey interactions consists of the following system of two differential equations:

\begin{matrix}\frac{dx}{dt} = P (x) − F (x, y); x(0) = x_{0} ≥ 0, \\ \frac{dy}{dt} = eF (x, y) − D(y); y(0) = y_{0} ≥ 0, \end{matrix}  (1.1)

where x and y represent the prey and predator population sizes, respectively, and functions P (x), D(y) describe the intrinsic growth rate of the prey and the mortality rate of the predator, respectively. The function D(y) is assumed to be linear (D(y) = dy), while the function P (x) has a linear (P (x) = ax) or logistic (P (x) = ax(1− \frac{x}{K})) expression. The function F (x, y) is called “functional response” or “feeding rate” and represents the prey consumption per unit time. The model assumes a linear correspondence between the prey consumption and the predator production through the positive constant e. Among the most popular functional responses used in the modelling of predator-prey systems are the Michaelis-Menten

2000 Mathematics Subject Classification. 37M05, 39A11, 65L12, 65L20. Key words and phrases. Finite-difference; nonstandard; positive; elementary stable; predator-prey; predator interference. \begin{matrix}c \\ ◯ \end{matrix}2007 Texas State University - San Marcos. Published February 28, 2007.

type F (x, y) = \frac{fxy}{x+c} and the ratio-dependent type F (x, y) = \frac{fxy}{x+by}. However, in some situations they predict unrealistic population dynamics of the predator or the prey. The Michaelis-Menten type does not account for the mutual competitions among predators [15], while the ratio-dependent type allows unrealistic positive growth rate of the predator at low densities [3, 13, 12]. The Beddington-DeAngelis functional response F (x, y) = \frac{fxy}{b+wx+y} was introduced independently by Beddington [2] and DeAngelis [6] as a solution of the observed problems in the classical predator-prey theory. It has an extra term in the denominator which models mutual interference between predators and avoids the “low densities problem” of the ratio-dependent type functional response.

The Beddington-DeAngelis predator-prey model with a linear intrinsic growth of the prey population, analyzed completely in [7], has the following nondimensional form:

\begin{matrix}\frac{dx}{dt} = x − \frac{axy}{1 + x + y}; x(0) = x_{0} ≥ 0, \\ \frac{dy}{dt} = \frac{exy}{1 + x + y} − dy; y(0) = y_{0} ≥ 0, \end{matrix}  (1.2)

where x and y represent the prey and predator population sizes, respectively, and the positive constants a, e and d represent the generalized feeding rate, generalized conversion efficiency and generalized mortality rate of the predator, respectively.

Numerical simulations of System (1.2) can be obtained by methods based on finite difference approximations, such as Euler, Runge-Kutta, or Adams methods. However, their use raises questions about the truncation errors, the stability regions and, from a dynamical point of view, the accuracy at which the dynamics of the continuous system are represented by the discrete system. The nonstandard finite difference schemes, developed by Mickens [16], have laid the foundation for designing methods that preserve the dynamical behavior of the approximated differential system. The techniques are based on a nonlocal numerical treatment of the right-hand side functions and more sophisticated discretizations of time derivatives. Nonstandard numerical techniques, preserving the stability of equilibria, have been applied to some specific dynamical systems by Lubuma and Roux [14] and Dimitrov and Kojouharov [9, 10], among others. The resulting so-called elementary stable nonstandard (ESN) methods solve the numerical stability problems locally for an arbitrary time step-sizes. However, the ESN methods, as well as the standard numerical methods, do not guarantee a positive discrete solution for all positive initial values. The positivity condition is natural when interspecies interactions are modeled and approximated numerically. Failing to satisfy it reflects negatively on the accuracy and efficiency of the numerical methods. Several classes of positive and elementary stable nonstandard (PESN) methods have been recently developed for some phytoplankton-nutrient [8] models and for a Rosenzweig-MacArthur predatorprey [11] model, in which a constant feeding rate per predator has been assumed.

The predator-dependent Beddington-DeAngelis functional response F (x, y) leads to richer dynamics of System (1.2), which makes the design and analysis of the PESN numerical methods more challenging that in earlier cases. The PESN methods, developed here, preserve both the positivity of the solutions and the stability of the equilibria of System (1.2). In addition, the designed numerical schemes allow for the discrete systems to be solved explicitly, which increases the efficiency of the methods.

The paper is organized as follows. Section 2 provides some definitions and preliminary results. Results from the mathematical analysis of the Beddington-DeAngelis predator-prey system followed by the development of the corresponding PESN numerical methods are presented in Section 3. In the last two sections we illustrate our theoretical results by a set of numerical examples and outline some future research directions.

2. Definitions and Preliminaries

A general n-dimensional autonomous system has the following form:

\frac{dx̄}{dt} = f(x̄); x̄(t_{0}) = x̄_{0},  (2.1)

where x̄ = [x^{1}, x^{2}, …, x^{n}]^{T} : [t_{0}, T) → ℝ^{n}, the function f = [f^{1}, f^{2}, …, f^{n}]^{T} : ℝ^{n} ↦ ℝ^{n} is differentiable and x̄_{0} ∈ ℝ^{n}. The equilibrium points of System (2.1) are defined as the solutions of f(x̄) = 0.

Definition 2.1. Let x̄^{∗} be an equilibrium of System (2.1), J(x̄^{∗}) be the Jacobian of System (2.1) at x̄^{∗} and σ(J(x̄^{∗})) denotes the spectrum of J(x̄^{∗}). An equilibrium x̄^{∗} of System (2.1) is called linearly stable if Re(λ) < 0 for all λ ∈ σ(J(x̄^{∗})) and linearly unstable if Re(λ) > 0 for at least one λ ∈ σ(J(x̄^{∗})).

A numerical scheme with a step size h, that approximates the solution x̄(t_{k}) of System (2.1) can be written in the form

D_{h}(x̄_{k}) = F_{h}(f; x̄_{k}),  (2.2)

where D_{h}(x̄_{k}) ≈ \frac{dx̄}{dt} |_{t=t_{k}}, x̄_{k} ≈ x̄(t_{k}), F_{h}(f; x̄_{k}) approximates the right-hand side of System (2.1) and t_{k} = t_{0} + kh. Throughout this article, we assume that System (2.1) has a finite number of hyperbolic equilibria, i.e., Re(λ) ≠ 0 for all λ ∈ Ω, where Ω = ∪_{x̄^{∗}∈Γ}σ(J(x̄^{∗})) and Γ represents the set of all equilibria of System (2.1).

Definition 2.2. Let x̄^{∗} be a fixed point of the scheme (2.2) and the equation of the perturbed solution x̄_{k} = x̄^{∗} + ϵ̄_{k}, where ϵ̄_{k} is small, be linearly approximated by

D_{h}ϵ̄_{k} = J_{h}ϵ̄_{k}.  (2.3)

Here the right-hand side of Equation (2.3) represents the linear term in ϵ̄_{k} of the Taylor expansion of F_{h}(f; x̄^{∗} + ϵ̄_{k}) around x̄^{∗}. The fixed point x̄^{∗} is called stable if ∥ϵ̄_{k}∥ → 0 as k → ∞, and unstable otherwise, where ϵ̄_{k} is the solution of Equation (2.3).

Let System (2.2) be expressed in the following explicit way:

x̄_{k+1} = G(x̄_{k}),  (2.4)

where the function G = [G^{1}, G^{2}, …, G^{n}]^{T} : ℝ^{n} ↦ ℝ^{n} is differentiable. If x̄^{∗} is a fixed point of System (2.4) then the equation for the perturbed solution ϵ̄_{k}, around x̄^{∗}, has the form:

ϵ̄_{k+1} = J(x̄^{∗})ϵ̄_{k},

where J(x̄^{∗}) denotes the Jacobian (\frac{∂G^{i}(x̄^{∗})}{∂x^{j}})_{1≤i,j≤n}. The solution ∥ϵ̄_{k}∥ → 0 when k → ∞ if and only if all eigenvalues of J(x̄^{∗}) are less than one in absolute values.

We introduce the next two definitions based on definitions given by Anguelov and Lubuma in [1].

Definition 2.3. The finite difference method (2.2) is called elementary stable, if, for any value of the step size h, its only fixed points x̄^{∗} are those of the differential system (2.1), the linear stability properties of each x̄^{∗} being the same for both the differential system and the discrete method.

Definition 2.4. The one-step method (2.2) is called a nonstandard finite-difference method if at least one of the following conditions is satisfied:

• D_{h}(x̄_{k}) = \frac{x̄_{k+1}−x̄_{k}}{φ(h)} , where φ(h) = h + O(h^{2}) is a nonnegative function;

• F_{h}(f; x̄_{k}) = g(x̄_{k}, x̄_{k+1}, h), where g(x̄_{k}, x̄_{k+1}, h) is a nonlocal approximation of the right-hand side of System (2.1).

The form of predator-prey systems (1.2) guarantees that any solution with positive initial conditions remains positive in time. The corresponding requirement for numerical methods is formulated in the following definition.

Definition 2.5. The finite difference method (2.2) is called positive, if, for any value of the step size h > 0, and x̄_{0} ∈ ℝ_{+}^{n} its solution remains positive, i.e., x̄_{k} ∈ ℝ_{+}^{n} for k = 1, 2, 3, …

The following theorem deals with the properties of the general nonstandard schemes based on standard finite-difference methods [1]:

Theorem 2.6. If the numerical method (2.2) represents a standard finite-difference scheme that is consistent and zero-stable, then any corresponding nonstandard finitedifference scheme in Definition 2.4 is necessarily consistent. Furthermore, if the nonstandard scheme is constructed according to the first bullet of Definition 2.4, then this scheme is zero-stable provided that the operator F_{h} satisfies, for some M > 0 independent of h and for any bounded sequences {x̄_{k}} and {x̃_{k}} in ℝ^{n}, the Lipschitz condition

sup_{k} ∥F_{h}(x̄_{k}) − F_{h}(x̃_{k})∥ ≤ M sup_{k} ∥x̄_{k} − x̃_{k}∥.

3. PESN Methods for Predator-Prey Models with Beddington-DeAngelis Functional Response

Depending on the values of the parameters in System (1.2) the following equilibria exist:

(1) E_{0} = (0, 0);

(2) E^{∗} = (x^{∗}, y^{∗}) = (\frac{ad}{ae−e−ad}, \frac{e}{ae−e−ad}). The equilibrium E^{∗} exists if and only

if ae − e − ad > 0.

According to the mathematical analysis of System (1.2) in [7], the following statements about the stability of the equilibria are true:

(1) The equilibrium E_{0} is always linearly unstable;

(2) The equilibrium E^{∗} is linearly stable if a < e and linearly unstable if a > e; The dynamics of System (1.2) are summarized in the following theorem [7].

Theorem 3.1. For the global dynamics of System (1.2) the following holds:

(1) If ae − e − ad ≤ 0 then all trajectories are unbounded. In addition, e > d yields that lim_{t→∞} x(t) = ∞ and lim_{t→∞} y(t) = ∞ while e ≤ d yields that lim_{t→∞} x(t) = ∞ and lim_{t→∞} y(t) = 0.

(2) If ae − e − ad > 0 and a < e then the interior equilibrium E^{∗} is globally asymptotically stable.

(3) If ae − e − ad > 0 and a > e then the interior equilibrium E^{∗} is unstable. If, in addition, e ≥ d + 1 then every trajectory circles indefinitely and counterclockwise around E^{∗} moving away from it. If e < d + 1 then for every trajectory is true that lim_{t→∞} x(t) = lim_{t→∞} y(t) = ∞.

(4) If ae−e−ad > 0 and a = e then the interior equilibrium E^{∗} is a global center, which means that all trajectories (except E^{∗} itself) are periodic orbits containing E^{∗} in their interior.

The new PESN methods for solving Beddington-DeAngelis predator-prey systems with linear intrinsic growth of the prey population can be designed, based on the following theorem:

Theorem 3.2. Let ϕ be a real-valued function on ℝ that satisfies the property:

ϕ(h) = h + O(h^{2}) and 0 < ϕ(h) < 1 for all h > 0.  (3.1)

Then the following scheme for solving System (1.2) represents a PESN method:

\begin{matrix}\frac{x_{k+1} − x_{k}}{φ(h)} = x_{k} − \frac{ax_{k+1}y_{k}}{1 + x_{k} + y_{k}}; \\ \frac{y_{k+1} − y_{k}}{φ(h)} = \frac{ex_{k}y_{k}}{1 + x_{k} + y_{k}} − dy_{k+1}; \end{matrix}  (3.2)

where φ(h) belongs to the following class of functions:

(1) If ae − e − ad ≤ 0 then φ(h) = ϕ(hq)/q for all q ≥ 0, with φ(h) = h for

q = 0.
(2) If ae − e − ad > 0 and a ≠ e then φ(h) = ϕ(hq)/q for q > \frac{e|a−2|}{|e−a|}.

Proof. The explicit expression of the nonstandard scheme (3.2) has the form:

\begin{matrix}x_{k+1} = \frac{(1 + φ(h))x_{k}}{1 + \frac{aφ(h)y_{k}}{1+x_{k}+y_{k}}} \\ y_{k+1} = \frac{(1 + \frac{eφ(h)x_{k}}{1+x_{k}+y_{k}} ) y_{k}}{1 + φ(h)d} . \end{matrix}  (3.3)

Since the constants a, e and d are positive then the scheme (3.3) is unconditionally positive and its fixed points are exactly the equilibria E_{0} and E^{∗} of System (1.2).

The analysis of the Jacobian J(E_{0}) of Scheme 3.3 at the equilibrium E_{0} shows that E_{0} is always an unstable fixed point of System (3.2).

The Jacobian of the Scheme (3.3) at the equilibrium E^{∗} has the following form:

J(E^{∗}) = (\begin{matrix}1 + \frac{φ(h)d}{(1+φ(h))e} & −\frac{φ(h)d(a−1)}{(φ(h)+1)e} \\ \frac{φ(h)(e−d)}{a(1+φ(h)d)} & 1 − \frac{φ(h)d}{a(1+φ(h)d)} \end{matrix})

The eigenvalues λ_{1} and λ_{2} of the Jacobian J(E^{∗}) are roots of the quadratic equation:

λ^{2} − (A + 2)λ + (A + B + 1) = 0,

where

A = \frac{φ(h)d}{(1 + φ(h))e} − \frac{φ(h)d}{a(1 + φ(h)d)}, B = \frac{φ(h)^{2}d(ae − ad − e)}{(1 + φ(h))(1 + φ(h)d)ae}.

For a quadratic equation of the form λ^{2} + αλ + β = 0 both roots satisfy |λ_{i}| < 1, i = 1, 2 if and only if the following three conditions are met: (a) 1 + α + β > 0; (b) 1 − α + β > 0; and (c) β < 1 [4, p.82]. Using this fact, we obtain that E^{∗} is a stable fixed point of Scheme (3.2) if and only if the following is true:

(a) B > 0;

(b) 4 + 2A + B > 0; and
(c) A + B < 0;

and that E^{∗} is an unstable fixed point if at least one of the above conditions fails.

The equilibrium E^{∗} exists when ae − ad − e > 0, which implies that a > 1 and e > d. Therefore B > 0 and the condition (a) is always true. Simple calculations yield that A > −1, which implies that the condition (b) is also true. The condition (c) is equivalent to the following inequality:

φ(h)e(a − 2) < e − a.  (3.4)

Assume that E^{∗} is a stable equilibrium of System (1.2). Therefore e − a > 0 and Inequality (3.4) is true for a ≤ 2. Let a > 2 and φ(h) = ϕ(hq)/q, where q > \frac{e|a−2|}{|e−a|}. Since 0 < φ(h) < \frac{1}{q} then φ(h) < \frac{e−a}{e(a−2)} for all h > 0. Therefore the condition (c) is satisfied and the interior equilibrium E^{∗} is a stable fixed point of the Scheme (3.2) with φ(h) = ϕ(hq)/q.

If the equilibrium E^{∗} is unstable then e − a < 0 and Inequality (3.4) is not true for a ≥ 2. This implies that the condition (c) fails and E^{∗} is an unstable fixed point of Scheme (3.2). If a < 2 and φ(h) = ϕ(hq)/q, where q > \frac{e|a−2|}{|e−a|}, then the condition (c) is not satisfied and the interior equilibrium E^{∗} is an unstable fixed point of the Scheme (3.2) with φ(h) = ϕ(hq)/q. □

The designed nonstandard methods (3.2) require that the dynamical system (1.2) has only hyperbolic equilibria and also that the parameter q is selected above a specific threshold value. In that case, the PESN methods replicate the properties (1), (2), and (3) of System (1.2) stated in Theorem 3.1, i.e., the methods guarantee that the stabilities of the equilibria E_{0} and E^{∗} of the differential system are the same as the stabilities of E_{0} and E^{∗} as fixed points of the numerical method for an arbitrary step-size h.

4. Numerical Simulations

To illustrate the performance of the designed new PESN finite-difference methods, we consider the Beddington-DeAngelis system (1.2) with a = 3.0, d = 2.25 and e = 4.0. Mathematical analysis of the predator-prey system shows that the equilibrium E_{0} = (0, 0) is unstable, while the equilibrium E^{∗} = (\frac{27}{5}, \frac{16}{5}) is globally asymptotically stable in the interior of the first quadrant [7].

The first set of simulations (Fig. 1 (a)-(c)) compares the performance of the PESN method (3.2), the explicit Euler method, and the ESN Euler method with θ = 0 [10]. For the above choice of system parameters the ESN method is elementary stable for q_{ESN} > 1.25 [10, Theorem 1] and the PESN method is elementary stable for q_{P ESN} > 4 (see Theorem 3.2), so we use q_{ESN} = 1.3 and q_{P ESN} = 4.5, respectively, in our set of numerical experiments. Numerical approximations of the solution of System (1.2) with initial conditions x_{0} = 5.0 and y_{0} = 3.5 show that for a relatively small step-size (Fig. 1 (a)) all three methods follow the correct dynamics of the original system, but an increase of the step-size (Fig. 1 (b)) causes the solution of the explicit Euler method, which is not elementary stable, to evolve away from the stable equilibrium E^{∗}. Moving the initial condition closer

\begin{matrix}(a) h = 0.1, x_{0} = 5.0, y_{0} = 3.5 & (b) h = 1.0, x_{0} = 5.0, y_{0} = 3.5 \\ (c) h = 1.0, x_{0} = 0.5, y_{0} = 6.5 & (d) h = 0.4, x_{0} = 7.5, y_{0} = 5.0 \\ (e) h = 1.27, x_{0} = 2.0, y_{0} = 3.5 & (f) h = 2.56, x_{0} = 2.0, y_{0} = 3.5 \end{matrix}

Figure 1. Numerical approximations of solutions to System (1.2).

to the x-axis, e.g. x_{0} = 0.5 and y_{0} = 6.5, causes the numerical solution of the non positivity-preserving ESN Euler method, after an initial jump of the solution’s trajectory outside of the first quadrant (Fig. 1 (c)), to evolve away from the equilibrium E^{∗}. The next simulation compares the approximations of the solution, originating at (x_{0}, y_{0}) = (7.5, 5.0), obtained by the PESN method (3.2) and the Patankar Euler method [5]. Even though the Patankar Euler method is positive, its numerical solution leads to a machine blowup for a variety of step-sizes (e.g. h = 0.4 in Fig. 1 (d)), because the method is not elementary stable. The last two simulations compare the PESN method with the second- and fourth-order Runge-Kutta methods. In the selected numerical examples the two Runge-Kutta methods, which are not elementary stable, fail to preserve the correct asymptotic behavior of System (1.2) for a variety of step-sizes (e.g. h = 1.27 and h = 2.56 in Fig. 1 (e) and (f), respectively).

Conclusions

The designed new PESN methods preserve two of the most important dynamical characteristics of predator-prey system with Beddington-DeAngelis functional response (1.2), namely the stability of all equilibria and the positivity of all solutions with positive initial conditions. In all of the numerical experiments (Fig. 1) the PESN solution converges toward the asymptotically stable interior equilibrium for arbitrary time-steps, while all of the other numerical methods fail to preserve either the stability of the equilibrium or the positivity of the trajectories for a variety of step-sizes. The solution of the PESN method follows the predicted by the mathematical analysis (Theorem 3.1) asymptotic behavior. It also preserves the local stability of the equilibria and the positivity of the solutions of System (1.2) for any step-size. The explicit form of the PESN methods makes them a computationally effective tool in the numerical simulations of predator-prey systems.

Future research directions include the design and analysis of PESN methods for the general predator-prey system (1.1) and the construction of similar nonstandard schemes for biological systems with non-hyperbolic equilibria.

References

[1] Anguelov R., Lubuma J. M.-S., Contributions to the mathematics of the nonstandard finite difference method and applications, Numer. Methods Partial Diff. Equations 17:5, (2001) 518-543.

[2] Beddington J.R., Mutual interference between parasites or predators and its effect on searching efficiency, J. Animal Ecol. 44, (1975) 331-340.

[3] Berezovskaya F., Karev G., Arditi R., Parametric analysis of the ratio-dependent predatorprey model, J. Math. Biol. 43, (2001) 221-246.

[4] Brauer F., Castillo-Chavez C., Mathematical Models in Population Biology and Epidemiology (Springer-Verlag, New York 2001)

[5] Burchard H., Deleersnijder E., Meister A., A high-order conservative Patankar-type discretization for stiff systems of production-destruction equations, Appl. Numer. Math. 47:1, (2003) 1-30.

[6] DeAngelis D. L., Goldstein R. A., O’Neill R. V., A model for trophic interaction, Ecology 56, (1975) 881-892.

[7] Dimitrov D. T., Kojouharov H. V., Complete mathematical analysis of predator-prey models with linear prey growth and Beddington-DeAngelis functional response, Appl. Math. Comput., 162:2, (2005) 523-538.

[8] Dimitrov D. T., Kojouharov H. V., Analysis and numerical simulation of phytoplanktonnutrient systems with nutrient loss, Math. Comput. Simulation, 70:1, (2005) 33-43.

[9] Dimitrov D. T., Kojouharov H. V., Nonstandard finite-difference schemes for general twodimensional autonomous dynamical systems, Appl. Math. Lett., 18:7, (2005) 769-774.

[10] Dimitrov, D. T. and Kojouharov, H. V., Stability-preserving finite-difference methods for general multi-dimensional autonomous dynamical systems, Int. J. Numer. Anal. Model., 4:2, (2007) 282-292.

[11] Dimitrov, D. T. and Kojouharov, H. V., Positive and Elementary Stable Nonstandard Numerical Methods with Applications to Predator-Prey Models, J. Comput. Appl. Math., 189:1-2, (2006) 98-108.

[12] Hsu S.-B., Hwang T.-W., Kuang Y., Global analysis of the Michaelis-Menten-type ratiodependent predator-prey system, J. Math. Biol. 42, (2001) 489-506.

[13] Kuang Y. and Beretta E., Global qualitative analysis of a ratio-dependent predator-prey system, J. Math. Biol. 36, (1998) 389-406.

[14] Lubuma J. M.-S. and Roux A., An improved theta-method for systems of ordinary differential equations, J. Differ. Equations Appl. 9:11, (2003) 1023-1035.

[15] May, R. M., Stability and Complexity in Model Ecosystem (Princeton Univ. Press, Princeton 1974)

[16] Mickens, R. E., Nonstandard finite difference model of differential equations (World Scientific, Singapore 1994)

Dobromir T. Dimitrov Department of Ecology and Evolutionary Biology, University of Tennessee at Knoxville, Knoxville, TN 37996-1610, USA

   E-mail address: ddimitr1@utk.edu

Hristo V. Kojouharov Department of Mathematics, University of Texas at Arlington, Arlington, TX 76019-0408, USA

   E-mail address: hristo@uta.edu