Open access peer-reviewed chapter

Numerical Simulation on Ecological Interactions in Time and Space

By Kornkanok Bunwong

Submitted: December 16th 2010Reviewed: May 9th 2011Published: September 9th 2011

DOI: 10.5772/24945

Downloaded: 2079

1. Introduction

The most formal and systematic tool to simplify a real-world phenomenon dealing with the interrelationships between organisms and their environment, including the interaction between each other is an ecological model. Of course, it is usually created for one or more purposes, for example, to gain system understanding, to forecast the future state of the system, and to develop new hypotheses. Not only that, a good model should also reach a balance between the complexities of the real-world system, which is too difficult to solve, and the need of simple formulation and valid analytic model, which can be explicitly solved. Ecological modeling has a long and rich history. So far many sophisticated ecological models have been developed in all fields of ecology such as the ecology of individuals including physiological ecology, the ecology of populations, and the study of ecosystems (classified by Hofbauer & Sigmund, 1988; Roughgarden, 1996). Ecological phenomenon can often be simplified by making some assumptions and studying it with suitable time scales and spatial interactions. It should be noted that a simple deterministic system could behave dramatically unlike when its time scale and its spatial interaction are changed. Therefore, the next coming sections aim to convince the significance of time scales and spatial interactions.


2. The classical concept of ecological modelling

The books on mathematical modeling in biology written by L. Edelstein-Keshet (1988), S. P. Ellner & J. Guckenheimer (2006), M. R. S. Kulenovic & O. Merino (2002), J. D. Murray (1993), and S. P. Otto & T. Day (2007) are the valuable teaching resources suitable for modelers who require theoretical concepts as well as desire mathematical and computational techniques to build and analyze models. In the early ecological models, their time scales were limited to either discrete or continuous. Usually, if populations of individuals have synchronized changes (i.e., reproduction, infection, recovery, migration, removal or mortality) at a regular time interval and no overlap between successive generations then the discrete time model is recommended. Otherwise, a continuous time model is preferred. Traditionally, a discrete time model is represented by using difference equation(s) whereas a continuous time model is constructed by using differential equation(s). Recently, the practical use of mathematical modeling is richly contributed. For example, L. J. S. Allen (1994), L. J. S. Allen & A. B. Burgin (2000), W. M. Getz & J. O. Lloyd-Smith (2006), and S. R.-J. Jang (2008) were interested in discrete epidemic models while K. Bunwong et al. (2009) studied nutrient removal process on a continuous time model.

According to spatial interactions, the early ecological models were based on the mass-action law, first coined in chemical reaction. Hence the system is homogenous and the spatial interaction neither exists nor plays an important role. Before the 1970s, mathematical modelers typically used ordinary differential equations, seeking equilibria and analyzing their stability (Neuhauser, 2001). Subsequently, exploring bifurcation diagrams and chaotic patterns has continuously received lots of attention, related articles were written by J. Awrejcewicz (1991), J. Awrejcewicz & C.-H. Lamarque (2003), and E. Ott (1993). In 1969, the implicit spatial model, known as metapopulation, was first introduced by R. Levins. He constructed a model to describe a population that consists of several sub-populations joined together with immigration and emigration. Levins's simple model was motivated by and applied to a pest control situation over a large region, within which local populations would fluctuate in asynchrony (Hanski & Gilpin, 1991). This new concept was closely linked with the processes of local extinctions and re-colonization. I. Hanski & M. Gilpin (1991) also provided a conceptual distinction between local, metapopulation, and geographical scales. Later on, I. Hanski (1999) applied this approach to conservation biology. However, the role of space at the individual level is still not directly mentioned.

Throughout this chapter, an SIS epidemic model, well known disease transmission model, is considered. Of course, there are lots of diseases in this world. Here we focus on the disease that does not produce immunity, for example, some STD’s, the eye disease, and the common cold. Therefore, the main situation is that the population is divided into a susceptible (S) group and an infectious (I) group. The S group is infected by the I group while the I group recovers from the disease and returns to the S group. Consequently, we set S(t)and I(t)to represent the number of S and I individuals at time t, respectively. The total population size is assumed to be constant, S(t)+I(t)=N. Moreover, (α/N)S(t)I(t)and γI(t)represent the infection rate at which the S population contracts the disease and the total number of I individuals who recover per unit time at the time t, respectively. Then, the SIS epidemic model should contain two equations. With constant population size, it can be reduced to a single equation. For continuous time scale, it is



In this case, a solution of Equation (1) always behaves non-oscillatory. For discrete time scale, the SIS epidemic model becomes


In this case, the behavior of the endemic solution can be very complicated as it can tend to an equilibrium point, to limit cycles, or show chaos. Obviously, time scale makes system behavior amazing. In the next section, an SIS model is extended with more different time scales. Finally, the qualitative structures are investigated.

3. The important of time scales

In more complex situation, flu virus can spread continuously while in season, mostly disappear for some period of time, and spread again in a new season. According to Fig. 1, the data shows periodic outbreaks of disease (the Influenza Division of Centers for Disease Control and Prevention, 2010). Obviously, the reasonable time scale for this model should be a combination of discrete and continuous scales. Thus the model either using difference equation(s) or differential equation(s) should be improper. Consequently, the calculus on time scales has been developed. S. Hilger first introduced this theory in order to unify continuous (R) and discrete (Z) analysis (Agarwal et al., 2002). Since then, the theory has been extended. Nowadays, time scales theory can be used to explain system dynamics not only for continuous and discrete times but also for other types of time such as period (P) and discrete jump with fixed length (hZ). There are plenty of publications in theoretical results. For example, J. Hoffacker & C. C. Tisdell (2005) studied on stability and instability for dynamic equations while E. Akin et al. (2001), D. R. Anderson (2009), and Y. Xu & Z. Xu (2009) studied on oscillation and nonoscillation criteria for dynamic equations and dynamic systems. However, there are few articles in numerical results (Sae-jie & Bunwong, 2009; Siming et al., 2008). In addition, there are some applications in economics (Atici et al., 2006) and epidemiology (Sae-jie et al., 2010a, 2010b; Thomas et al., 2009). In order to contribute time scales concepts and techniques understandably, some useful notations and definition are introduced as follows.

Figure 1.

The number of influenza-associated pediatric deaths.

3.1. Notations and definitions

A time scale, T, is an arbitrary nonempty closed subset of the real numbers. Forward and backward jump operators are defined by σ(t)=inf{sT:s>t}and ρ(t)=sup{sT:s<t}, respectively where inf=supT, sup=infT, and denotes the empty set. A point tTis called left-dense if t>infTand ρ(t)=t, right-dense if t<supTand σ(t)=t, left-scattered if ρ(t)<t, right-scattered if σ(t)>t, isolated if ρ(t)<t<σ(t), and dense if ρ(t)=t=σ(t). Moreover, the set Tκis defined by T\{m}if Thas a left-scattered maximum m. Otherwise, it is T. Finally, the graininess function μ:T[0,)is defined by μ(t)=σ(t)t(Bohner & Peterson, 2001).

3.2. Calculus on time scales

Traditionally, limits and continuity are key concepts for calculus development including calculus on time scales. A function f:TRis said to be rd-continuous (right dense continuous) provided fis continuous at right-dense points and left-hand limits exist and it is finite at left-dense points in T. Assume f:TRis a function and tTκ. Then the following statements are equivalent:

The (delta) derivative of f:TRat point tTκ, fΔ(t), exists.

For all ε>0, there is a neighborhood Uof t(i.e., U=(tδ,t+δ)Tfor some δ>0) such that for all sU,


Apart from discrete and continuous time scales, there are two more interesting time domain (as visualized in Fig. 2) to introduce. The first one is the combination of continuous and discrete time scales, denoted by the symbol T=Pl,h:=k=0[k(l+h),k(l+h)+l]where l,h>0and kN0. For this period time scale, lis the fixed length of the continuous interval while his the fixed length of the discrete jump. The second one is composed of points that are equally spaced in time. Suppose that the distant between two successive points is h. Therefore, the symbol of this time scale is T=hZ={hk:kZ,h>0}. Later on, forward and backward jump operators, graininess function, and (delta) derivative of function for each time scale are provided.

Figure 2.

Examples of time scales.

In the case T=R, we have

σ(t)=ρ(t)=t, μ(t)=0,


In the case T=hZ, we have

σ(t)=t+h,ρ(t)=th, μ(t)=h,


In the case T=Pl,h, if L=k=0[k(l+h),k(l+h)+l), H0=k=0{k(l+h)+l}, and H1=k=1{k(l+h)}then we have





For more details on time scales, we refer the reader to Bohner & Peterson (2001,2003)

3.3. SIS epidemic model on time scales

After replacing dSdtin Equation (3) with SΔ(t), we obtain the single equation of SIS epidemic model on time scales. Previously, one solution of continuous time model tends to one asymptotically stable equilibrium point while the solution of discrete time model displays various behaviors depending on parameter values. The following analyzes reveal system behavior on different time scale via two approaches. Firstly, we are interested in SIS epidemic model on T=hZ, i.e.,


It should be pointed out that μand hare equivalent from now on. μis not a graininess function. Obviously, T=Zwhen h=1and T=Rwhen htends to zero. Secondly, we change our focus to SIS epidemic model on T=Pl,h, i.e.,




Then we explore the numerical results of previous dynamic equations individually. For T=hZtime scale, we particularly investigate the system behavior when hor μvary. For T=Pl,htime scale, we fix h=1. Therefore, T=Zwhen lapproaches zero.

3.4. Numerical solution

Since the numerical method for ordinary differential equation and difference equation are already well-known, this section contains only some additional algorithm for computing the numerical solution of dynamic system on four time scales (Sae-jie, & Bunwong, 2009). Before running the following program, the user must define the function F and G.

INPUT: A: starting time

X, Y: initial values

M: number of step size (Case 2: M = l/H)

H: step size for continuous interval

Case 1:


FOR j = 0, 1, 2,..., M DO

T = A + H*j




Case 2:


P: numbers of period

l: length of continuous interval

h: length of jump

DUM3 = A

FOR j = 0, 1, 2,..., M DO

T = DUM3 + H*j


DUM1 = X

DUM2 = Y



X = DUM1 + h*F(T, DUM1, DUM2)

Y = DUM2 + h*G(T, DUM1, DUM2)

FOR k = 1, 2, 3,..., P DO

DUM3 = A + k*(l+h)

FOR j = 0, 1, 2,..., M DO

T = DUM3 + H*j


DUM1 = X

DUM2 = Y



X = DUM1 + h*F(T, DUM1, DUM2)

Y = DUM2 + h*G(T, DUM1, DUM2)


Case 3:


DUM3 = A

FOR j = 0, 1, 2,..., M DO

T = DUM3 + H*j


DUM1 = X

DUM2 = Y

X = DUM1 + H*F(T, DUM1, DUM2)

Y = DUM2 + H*G(T, DUM1, DUM2)


F := XΔ(T, X, Y)

G := YΔ(T, X, Y)

In the case T=hZ, W. Sae-jie et al. (2010a) always fixed the following parameter values γ=0.9and N=100. For each value of α, μwas treated as the bifurcation parameter. When α=2, the dynamic behavior of continuous and discrete system are the same. For α=3.4, the solutions, however, appear as asymptotically stable, a period two cycle, a period four cycle when μ=0.1,0.9,1.0, respectively as shown in Fig. 3.

Figure 3.

The time series solution ofEquation (9)whenα=3.4.

M. R. S. Kulenovic & O. Merino (2002) and L. Tien-yien & J. A. Yorke (1975) discovered that if there exists a period three cycle, then there exists chaotic behavior. For α=3.6, the bifurcation diagram (as shown in Fig. 4) exhibits a period three cycle for μ(1.0476,1.0524). Consequently, dynamic equation (9) can generate chaotic pattern for a proper value of μ. For more details, the non-oscillatory solution, the oscillating period two solution, and the chaos occur when μ=0.1, μ=0.8, and μ=1, respectively. The last behavior is illustrated in Fig. 5.

In the case T=Pl,h, this period time scale and a continuous time scale are equivalent when μ=h=0. Thus, the result appears as a non-oscillatory solution. However, W. Sae-jie et al. (2010b) always fixed the length of discrete jump together with the following parameter values α=3.6, μ=h=1, γ=0.9, and N=100but varied the length of continuous interval, l. For sufficiently high value of l, the result (as shown in Fig. 6) appears as a non-oscillatory solution which is similar to the result in Fig. 3. for a continuous time scale. It disappears in some intervals because of the discrete time jump. When the length of the continuous interval decreases, a period two cycle with some continuous intervals appears as visualized in Fig. 7. If the length of the continuous interval is gradually reduced and closely zero, T=Pl,his similar to a discrete time scale.

Figure 4.

The bifurcation diagram ofμwhenα=3.6.

Figure 5.

The time series solution ofEquation (9)whenα=3.6andμ=1.

4. The important of spatial interaction

In contrast to a metapopulation model, space in an explicit model should be taken into account as another variable. As for non spatial models, the time scale, space, and population state for an explicit model can be either discrete or continuous. If all of them are continuous variables, a partial differential equation has been used (for example, Fisher, 1937). However, solving partial differential equations analytically is difficult. Consequently, analysis is based on computer simulation. Then we face another limitation because solving them numerically

Figure 6.

The time series solution ofEquation (10)on time scaleP1,1.

Figure 7.

The time series solution ofEquation (10)on time scaleP0.1,1.

involves discretizing space and time. In principle, this reduces them to coupled map lattice model. On the other hand, if space is treated as discrete, a lattice model is proposed. With respect to state variables, lattice models can be divided into three subgroups, these are coupled map lattice (Morris, 1997), cellular automaton (Wolfram, 1986), and network model (Sole & Manrubia, 1996, 1997; Verdasca et al., 2005). A major problem with the lattice model is that analysis is often restricted to direct computer simulation which consumes expensive time when it is a stochastic model on a reasonably large lattice.

Fortunately, in microscopic point of view, there is an alternative way, so-called pair approximation, that can help us not only to produce the numerical solutions of spatially explicit lattice models by using ordinarily numerical method but also to analyze them. The Japanese researchers were pioneers to apply this idea to ecological systems and have continuously produced theoretical result and applications (Harada et al., 1995; Iwasa et al., 1998; Kubo et al., 1996; Matsuda et al., 1992; Sato et al., 1994). They referred to the pair approximation as the doublet decoupling approximation and represent their system as coupled system of ordinary differential equation for the density of singletons and pairs. However, Rand’s colleagues (Bauch, 2000; Bunwong 2006; Keeling, 1995; Morris, 1997; van Baanlen & Rand, 1998) used pair approximation as a moment closure approximation where the number of pairs and the number of singletons are the only state variables. Higher order correlations are neglected. A more advantage of this method is that it keeps spatial correlations. There are varieties of applications in ecological interactions (Dieckmann et al., 2000; Ellner, 2001; Ellner et al., 1998), epidemiology (Benoit et al., 2006; de Aguiar et al., 2004 ; Dieckmann et al., 2002 ; Elliott et al., 2000 ; Joo & Lebowitz, 2004), and forest dynamics (Schlicht & Iwasa, 2007). However, the mathematical formulas are still limited. Thus, we attempt to broaden the pair approximation idea by calculating configuration averages. We call our method as a new approach and previous method as an original approach. In order to contribute our techniques understandably, some useful notations and definition are introduced as follows.

4.1. Notations and definitions

Under a given configuration σ¯=(σk)where k{x,e}, the following Rand’s notations are defined (Rand, 1998).

σx,σeare the state of the site xand the edge e, respectively,σx=iis that the state of the individual xis i, σe=ijis that one end of the edge eis in state i, ei, while the other is in state j, ej, [i],[ij],[ijk]are the number of sites, edges, and triples in state i, ij, and ijk, respectively, Qx(i),Qej(i)are the number of i-state neighbors of the sites xand ei, respectively,Qx(i)σx=jis the average value of the number of i-state neighbors of a j-state site, Qej(i)σe=jkis the average value of the number of i-state neighbors of a j-state site in a jk-state edge,qiequals [i]/Nwhere Nis the total population size, andqi|jequals [ij]/Q[j]where Qis the average number of neighbors.

In this framework, space is represented by a network of sites. Each site can either be occupied by an individual or remains as an empty site that is still available for an individual to occupy. Two sites are neighbors when they regularly interact with each other. Joining these two neighboring sites performs an edge or pair. A line is used for this interaction. Fig. 8 provides an example when the state of site xis iand the state of site yis j. Both sites are neighbors. Moreover, there is an edge e.

Figure 8.

4.2. SIS Master and correlation equations

Our correlations are microcorrelations which can be measured on the scale of the interactions of individuals. After approximating higher order terms in master equations, we obtain a system of ordinary differential equation which composed of density of lower order terms, known as correlation equations. Moreover, the approximation technique is called the moment closure approximation. The differential equation for the single numbers involves pair numbers, triple numbers, etc. The differential equation for the pair numbers involves triple numbers, etc. So we get an infinite hierarchy of equations. However, we need to truncate the hierarchy at some point. For instant, pair approximation, the first order of moment closure approximation, truncates triples and higher order terms as functions of singletons and pairs only (Rand, 1998).

Let fbe a real-valued function of the state of the network at time t, which can be approximated as continuous. The equation fis derived by summing over all events in the population which affect fand the total change produced by those events is


where dfdt=εeventsr(ε)Δfεis the rate of event r(ε)and εis the change produced in fby event Δfε. It is called the master equation.

For our case study, the state of each site and edge will change over time as a consequence of two major types of events – infection and recovery. Infection changes the state εof the edge σe=SIinto the state eat rate σ'e=IIand recovery changes the state βof a site σx=Iinto the state xat rate σ'x=S. Therefore, the SIS spatial model becomes


If d[II]dt=2σe=SIβQeS(I)2σx=IδQx(I)and δare constant, the original approach of pair approximation is still valid. However, the real-world situation is more complicated. For example, the human-to-human transmission of Swine Flu occurs by inhalation of infectious droplets and droplet nuclei, and by direct contact, which is facilitated by air and land travel and social gatherings (Sinha, 2009). Therefore, the transmission rate and the recovery rate could vary depending on the surrounding infectious people. Consequently, we are able to assume that the infection rate and the recovery rate are βand β=b0+b1QeS(I), respectively where δ=d0d1Qx(I)are constant. Then the formulation of pair approximation is in trouble. Using the fact that



we, then, can use average forms instead of summation terms in the master equation as follows,


4.3. Calculating configuration averages

It should be pointed out that there are two types of average values in this new approach. Previously, all mentioned average values are called space or population average value because they are averages of the quantity over subsets of population. Now we introduce another way to calculate an average value. That is the expected value with respect to probability distribution. This seems reasonable to assume that if population size is large enough, then the configuration averages approximate probability expectations. Under some necessary and sufficient conditions, the space average and probability average are identical. K. Bunwong (2006, 2010a, 2010b) developed more formulas under multinomial distribution with parameters d[II]dt=2b0[SI]QeS(I)σe=SI+2b1[SI]QeS(I)QeS(I)σe=SI2d0[I]Qx(I)σx=I+2d1[I]Qx(I)Qx(I)σx=Iand Qwhere piand Poisson distribution for coordination numbers, respectively. In this chapter, the SIS spatial model is based on the framework that each site connects to a fixed number of neighbors, pi=qi|j. The following Bunwong’s formulas (2006) are used.


4.4. Numerical solution

The numerical method for spatial model (14) is as same as for ordinary differential equation. Here, the number of pairs and the number of singletons are state variables. K. Bunwong (2010a) mainly investigated the density of infected individuals, defined by Qej(l1)Qej(l2)σe=ij={(Q1)ql1|j+(Q1)!(Q3)!ql1|j2;l1=l2(Q1)!(Q3)!ql1|jql2|j;l1l2=qi, along the time series. We always fixed [I]/N. In case that the infection rate and the recovery rate are not affected by the surrounding infectious individuals (b0=0.3,d0=0.2), Fig. 9 reveals that the more nearby neighbors, the higher density of infected individuals at the equilibrium point. Moreover, in case that only the infection rate is affected by the surrounding infectious individuals (b1=0,d1=0), the stronger effect of the surrounding infectious individuals on the infection rate, the higher density of infected individuals at the equilibrium point. In case that only the recovery rate is affected by the surrounding infectious individuals (b10,d1=0), the stronger effect of the surrounding infectious individuals on the recovery rate, the higher density of infected individuals at the equilibrium point. The illustrations can be seen in Bunwong (2010a).

Figure 9.

Time evolution of the density of infected individualsb1=0,d10. Parameters:qiandb0=0.3,b1=0,d0=0.2,d1=0(from top to bottom, respectively).

5. Conclusion

Obviously, the assumption on timescale affects the system behavior. In order to simplify real-world correctly, choosing a suitable time scale is important for model formulation. Our example, SIS epidemic model, proves that if continuous time scale is used then two solutions of the system are asymptotically stable or unstable depending on parameter values and stable oscillating solutions have never existed. In contrast, if discrete time scale is applied then there are various types of solution behaviors such as equilibrium point solutions, period two cycles, period four cycles, period three cycles, and also chaotic solutions depending on parameter values as well. Consequently, the predicted behaviors from a model can be qualitatively very different for different time scales. The hard to answer question is “What is the proper model to understand observed data in a variety of time measurements?”. Of course, the observed data is usually discrete. Can we use differential equation(s)? With theory of time scales, the models are more varieties. Bifurcation diagram shows that the discrete jump and the continuous interval are essential. The differential equation(s) can be described observed data if the jump distant of data is sufficiently small. The distant also depends on other parameters. Moreover, numerical verification can visualize other interesting behavior pattern and guide other mathematician to consider theoretically.

In recent years, the effect of spatial structure is often taken into account in ecological interactions. Pair approximation is one of the powerful tools to understand human to human interactions. We have developed the way to approximate higher order quantities and applied to ecological problems. Particularly, our new approach is suitable for a model evolving according to the transition rates affecting additionally by neighbors. For example, people infect flu virus easily from their nearby neighbors. The health organization usually suggests infectious people to have some rest and be away from public places. It implies that if we surrounding with more infectious people, then we have higher chance to infected and/or lower chance to recover as shown in the numerical results.

Further studies on time scales should involve more numerical techniques and applications while future works on spatial interaction should concern development of formulation.



I would like to thank Prof. David Rand for introducing me to spatial world, Prof. Yoh Iwasa for sharing his experience, Prof. Ravi Agarwal for introducing me to time scale world, Dr. Elvin Moore for intuitive ideas, and Dr. Wichuta Sae-jie for collaboration. I am grateful to the Development and Promotion of Science and Technology Talent Project (DPST) and Thailand Research Fund (TRF) for financial supports under contact number MRG5180246.

© 2011 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike-3.0 License, which permits use, distribution and reproduction for non-commercial purposes, provided the original is properly cited and derivative works building on this content are distributed under the same license.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Kornkanok Bunwong (September 9th 2011). Numerical Simulation on Ecological Interactions in Time and Space, Numerical Analysis - Theory and Application, Jan Awrejcewicz, IntechOpen, DOI: 10.5772/24945. Available from:

chapter statistics

2079total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Unscented Filtering Algorithm for Discrete-Time Systems with Uncertain Observations and State-Dependent Noise

By R. Caballero-Águila, A. Hermoso-Carazo and J. Linares-Pérez

Related Book

First chapter

Application of the Lyapunov Exponents and Wavelets to Study and Control of Plates and Shells

By J. Awrejcewicz, V.А. Krysko, I.V. Papkova, Т.V. Yakovleva, N.A. Zagniboroda, М.V. Zhigalov, A.V. Krysko, V. Dobriyan, E.Yu. Krylovа and S.A. Mitskevich

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us