Cellular Automata in Ecological and Ecohydraulics Modelling

P.M.A. Sloot, B. Chopard, and A.G. Hoekstra (Eds.): ACRI 2004, LNCS 3305, pp. 502–512, 2004.© Springer-Verlag Berlin Heidelberg 2004Cellular Automata inEcological and Ecohydraulics Modelling Arthur Mynett and Qiuwen ChenWL | Delft Hydraulics, P.O. Box 1772600 MH Delft, The Netherlands{arthur.mynett, qiuwen.chen}@wldelft.nlAbstract. Cellular Automata are discrete dynamical systems in which many simple components act together locally to produce complex patterns on a global scale, which may exhibit “self-organising” behaviour. Owing to the ability to model local interactions and spatial heterogeneity, cellular automata have been applied to very broad fields. This paper presents the application of cellular automata to modelling (i) a confined ecosystem, (ii) prey-predator population dynamics, where evolution rules are defined entirely by geometric relations,and (iii) open aquatic ecosystems where external forcings are accounted for in the definition of cell state transitions. The results indicate that cellular automata could be a valuable paradigm in ecological and ecohydraulics modelling.1 IntroductionCellular automata (CA) constitute a mathematical system in which many simple components act together to produce complicated patterns of behaviour. CA models define interactions between species at a local level. Complex interactions emerge on a global scale through the evolution of relatively simple local rules. Such global behaviour may be very difficult – perhaps even impossible – to describe by a set of deterministic equations due to complicated temporal and spatial relationships. The ‘self-organising’ property of CA implies that a randomly distributed population can generate ordered structures through a process of irreversible evolution (Chen, 2004).Cellular automata have been applied in very broad fields, such as lattice gas hydrodynamics, urban dynamics and forest fire spreading. Owing to the ability to deal with local interactions and spatial heterogeneity, cellular automata may provide a viable approach in ecological modelling. This paper presents the capabilities of cellular automata for modelling ecological and ecohydraulics systems through several case studies. Applications include macrophyte vegetation dynamics (Chen, et al.,2002), population dynamics (Chen and Mynett, 2003) and harmful algal bloom prediction (Chen and Mynett, 2004). The model results indicate that CA could be a valuable alternative for differential equations in ecological / ecohydraulics modelling.2 Cellular Automata for Modelling Prey Predator DynamicsModelling population dynamics is of increasing interest in ecological research. These kinds of ecological models are used either for qualitative estimates about species extinction or for more quantitative estimates of spreading and interaction of multiple 网络自动机模拟不连续的，离散的动力系统展现出自组织由于交互作用空间异质性把…应用到封闭的生态系统捕食进化规则，进化论清晰的几何关系开放的水生系统外力证明定义细胞状态转换显示有价值的范例建立定义交互作用出现全球尺度相对简单地确定性方程复杂的时间，空间关系性质随机分布的人口不可逆演化产生有序结构应用到晶格气流体动力学城市动态森林火灾蔓延可行的方法水生植物赤潮选择微分方程定性评估物种灭绝定量评估Cellular Automata in Ecological and Ecohydraulics Modelling 503species under various environmental conditions. They also provide decision support to natural resources management like sustainable harvesting strategies in fishery industries. Population dynamics, however, is not very easy to describe mathematically due to the high complexity of the spatial and temporal relationships. Conventional population dynamics models are mostly aggregated, lumping species into biomass and formulating the dynamics by means of ordinary or partial differential equations. Oneof the simplest paradigms is provided by the predator-prey equations of Lotka-Volterra (LV). However, these models may fail to produce realistic results when individual properties and local interactions play a significant role in determining the relationships between populations and between species and its surroundings.An alternative approach is to apply individual-based modelling techniques or spatially explicit modelling paradigms. In this paper, a cellular automata based prey-predator model (EcoCA) is developed and compared with the Lotka-Volterra model.2.1 Development of EcoCAEcoCA is a stochastic cellular automata based computer simulation model that was developed to approximate a simple prey-predator system (Chen, 2004). It is a two-dimensional CA model which has up to three possible cell states: empty, prey and predator. The boundary conditions are fixed in such a way that the neighbourhoods are completed with cells taking the state of empty. The state of each cell is exclusive;namely that at each time step only one of the three states can exist in one cell. Theevolutions for each cell (i, j) are based on cell current state S i j t ,, the number of itsneighbouring cells occupied by prey, N p yt, and the number of its neighbouring cellsthat are occupied by predator, Npdt , as given by,(,,)t t ti j py pd p f S N N =(1)where f are evolution rules defining the transition probability p that a cell will becomeeither prey or predator dominated or empty at the next time step. These evolution rules take into account reproduction, mortality, predation, loneliness and overcrowdedness. Snapshots of a random initial condition and the spatial pattern at t =600 of a demonstrated run are presented in Fig 1.Fig. 1. Snapshots of initial condition (t = 0) and population dynamics (t= 600)不同环境状况为支持天然资源管理可持续开采策略海捞法规人口动力空间时间传统的大多聚集聚集表达普通或部分微分方程捕食方程？现实的个人财产替代方法个人本位随机的模拟模型接近二维空间的最多具有单元状态边界条件邻域唯一的换句话说跃迁概率，转变几率支配复制死亡捕食行为？空满快照随机初始条件空间图形504 A. Mynett and Q. ChenSeveral runs of this system with various starting conditions all lead to very similar patterns and results as exemplified by the plot of the population dynamics and phase dynamics (Fig 2).Fig. 2. Population dynamics of prey and predator (left), and phase trajectory (right)Table 1. Statistical descriptors of EcoCA model outputMean Standard deviationSkewness Prey 38972.20.11Predator10618.40.112.2 Quantification of Model ResultsThe model is tested to be trend-free, and it also passes the F -test that the variance is constant as well as the t -test that the mean is constant, by using split-record tests. This means that the system is stationary and the simple statistical descriptors given in Table 1 can describe the outputs from the CA ing this way to quantify and visualise the performance of the CA model, it becomes possible to compare the CA model results with real and measured data. As a starting point it was decided to produce a test data set based upon the classical spatially homogeneous LV prey-predator model. This simple deterministic model provides a continuous description of a two-species competition system and may be described in terms of the following ordinary differential equation system:2dGaG bG GH dt α=−−(2)dHcH GH dt β=−+(3)where G represents the population of prey, H represents the population of predators, a is the intrinsic growth rate of prey, b is the carrying capacity limit, c is the mortality rate of predator, α is the functional response of prey on predator, and β is the functional response of predators to prey. Each and every organism in the prey population (G ) and in the predator population (H ) is considered only in relation to thetotal quantity of its functional type. The effect of the predators upon the prey population is only measured by the ‘functional response’ term GH .可作榜样的情节相轨迹量化无趋势变化恒定的拆分记录不动的说明输出的信息形象化同质的捕食模型连续的形容以…的观点常微分方程被掠食者的数量捕食者的数量被掠食者本身的增长速率承载能力极限捕食者的死亡率功能反应生物，有机体与…有关功能型种群动力Cellular Automata in Ecological and Ecohydraulics Modelling 505These equations are highly non-linear. A 4th order Runge-Kutta method was applied to solve these equations numerically, using a time step of 0.2 days. In order to examine whether the CA model output was at all representative for the population dynamics described by the LV equations, the LV model was first calibrated to minimise the mean square error between the CA model output and the LV model output. The population dynamics and phase trajectory of the LV model is given in Fig 3.Fig. 3. Population dynamics of prey and predator (left), and phase trajectory (right)The LV model obviously displays a more uniform and regular periodicity. This is of course a result of the purely deterministic nature of the governing equations, such that the period of oscillation is largely determined by the parameters of the model and the amplitude is determined solely by the initial conditions. Although severe reservations can be expressed about the realism of the LV model for modelling natural systems, the resulting stable limit cycle is an observable property of real populations (May, 1976). Nevertheless, in the absence of real data, the results for theLV model do provide us with a set of data that can be used to perform a quantitative comparison with the CA model. Using the methodologies described above we can quantify the magnitudes of the trend and the periodic components. Once more, the time series was tested for the absence of trend and the invariability of the mean and variance. Having satisfied these tests, the statistical properties of the LV time series were then calculated. These properties are summarised in Table 2.Table 2. Statistical descriptors of LV model outputMean Standard deviation Skewness Prey 389550Predator112140Comparing the values of the statistical descriptors in Table 2 with those in Table 1, it is seen that the means of the two populations are very similar in both cases. The variance of the CA model, however, is significantly larger than the LV model and, of course, the LV model results are completely uniform and thus have a skewness of zero. The differences in the CA model can be explained by the stochastic componentin the CA model, which is entirely absent in the LV model.Based on the results of this非线性的适用于不惜一切？校准求最小值平均平方误差种群动力相轨迹统一的周期性规律的基本方程振动振幅单独地初始条件真实性由此产生的稳定极限环可观测的财富然而缺席比较趋势的幅度周期分量不变性均值和方差概括均值标准偏差偏斜度数据说明方差完全均匀随机分量506 A. Mynett and Q. Chenanalysis, it seems fair to say that the mean field results of the CA model quite closelymatch the population dynamics represented by the classical LV model output.3 Cellular Automata for Modelling Macrophyte Dynamics Clearly, the dynamics of Conway’s “Game of Life” and the EcoCA model introduced in section 2 is still purely dependent on geometric relations between neighbouring cells, and does not account for physical or biological processes. Therefore, the approach followed here is to extend the conventional cellular automata paradigm by incorporating external forcing. The CA model dynamics are thus not solely governed by the geometric relations of the particular stencil that was chosen, but also by external factors (Chen, et al., 2002a), as illustrated in Fig 4. The first case study involves a CA based model developed to simulate the competition and succession oftwo macrophyte species in the eutrophic Lake Veluwe, the Netherlands.Fig. 4. Diagram of CA model coupling with other models3.1 Description of Study AreaLake Veluwe is an artificial isolated part of the larger Lake IJssel in the centre of the Netherlands. The total water surface is around 3300 ha, with an averaged depth of 1.4m. It was formed by the construction of dams in the Southeast part of Lake IJssel in 1952 (Fig 5). According to long-term documentation, the submerged vegetation ofthe lake has experienced a great change after its formation due to the change in nutrient loading. Before 1968, the water in the lake was clear, with diverse macrophyte vegetation. Due to discharge of wastewater from some small cities, the lake was eutrophicated, and blue-green algae became dominant. Some restoration measures were taken in late 1970s, which resulted in the increase of P. pectinatus .The increase of P. pectinatus provided the precondition for the return of C. aspera .After 1990, C. aspera colonised steadily and gradually replaced the dominance of P.pectinatus .From an ecological point of view, it seemed that P. pectinatus would outcompete C. aspera in this lake system. However, C. aspera outcompeted P. pectinatus and replaced it gradually in Lake Veluwe. Analysis of long-term observations indicated a self-reinforcing ability of C. aspera during eutrophication. C. aspera returned at a lower phosphorus level (0.1 mg/l) than the level at the time of their disappearance (0.3平均场结果水生植物动态几何关系对…负有责任延长，扩展传统的网格自动机范例结合外部力单独地几何关系模板模仿人工的孤立的公顷通过…形成水坝建设水中浸没的植被形成营养负荷不同的水生植被污水排放富营养化蓝绿藻恢复措施先决条件稳步拓殖在与…竞争中胜出长期观察指出自增强能力富营养化磷水平Cellular Automata in Ecological and Ecohydraulics Modelling 507mg/l), a phenomenon known as hysteresis, therefore phosphorus is not a key factor in this case. It is supposed that the competition of dissolved inorganic carbon -13HCO and competition of light are the two main factors of the succession. However, the replacement process remained unclear, and triggered the demand for model simulation. In view of the environmental heterogeneity and the local interactions, this research selected a CA approach to simulate the competition of light and -13HCO in order to explain the essential features of the replacement process.Fig. 5. The study area Lake Veluwe3.2 Model DevelopmentIn this CA model, deterministic evolution rules obtained through laboratory and field experiments are applied. The model is designed to contain two partly interacting parallel submodels, one for P. pectinatus and the other for C. aspera . The processes considered in each submodel include shading, attenuation, -13HCO competition,photosynthesis, respiration, morality and spreading. A conceptual framework of the model can be found in Chen, et al., (2002). General aspects of the model include: (1)germination of P. pectinatus and C. aspera from propagules; (2) initialisation with exponential growth rate; (3) growing and spreading; (4) production of propagules.The deterministic evolution rules used in this CA model are based on laboratory experiments and are calibrated against field observations. The algorithms involved are presented in some detail in Chen et al. (2002).3.3 Results and DiscussionThe simulations of the model are presented in two ways: by visualisation of the growing and spreading patterns of the two species in the lake, and by time series of biomass density averaged over sampled cells. Two graphs (Fig 6) show the changes of biomass density of P. pectinatus and C. aspera in each cell and the colonisation process resp. As shown in Fig 6, the colonisation is from the Northeast to the Southwest, and it is faster in longitudinal direction than in transverse direction.Besides, the colonisation of C. aspera is faster than that of P. pectinatus . SeveralVeluwe滞后现象磷溶解无机碳猜想引发对模型仿真的需求异质性通过实验室获得适用交互并行子模型光合作用呼吸作用概念框架一般特性发展繁殖体初始化指数增长速率用…进行标定运算法则可视化生物量密度样品池定植纵向横向508 A. Mynett and Q. ChenC. asperaP. pectinatusC. asperaP. pectinatussimulation scenarios were carried out to test the governing factors. The results showedthat the light intensity and -13HCO are two major factors to the competitive growths of C. aspera and P. pectinatus in Lake Veluwe. Thus, shading and competition of-13HCO become two important processes. The scarcity of -13HCO has a great negative effect on the growth of P. pectinatus , while it has a indirect positive effect onthe growth of C. aspera , which is advantageous to the replacement of P. pectinatus by C. aspera . These results are compatible with the field observations of Marcel (1999),who explored an individual based model to study the effects of changing environmental conditions on the dynamics between C. aspera and P. pectinatus . It can be concluded that a CA model with deterministic local interaction rules proved to be a good modelling paradigm where other methods had failed so far.Fig. 6. Germination of both species (P. pectinatus germinates earlier than C. aspera ) at the left,and colonisation pattern by the end of the second year (P. pectinatus spreads slower than C.aspera ) at the right4 Fuzzy Cellular Automata for Modelling Harmful Algal BloomsThe EcoCA and the macrophyte models use either stochastic rules or deterministic rules. However, there are many systems, in particular ecosystems such as harmful algal blooms, where the detailed mechanisms and their statistical behaviour remain unclear. For that reason, a fuzzy logic technique is introduced into the cellular automata paradigm for rule formulation, denoted as Fuzzy Cellular Automata.4.1 Description of the Study Area The Dutch coast receives outflow from the rivers Rhine and Meuse and is one of the most productive fishing areas in the world. In the past 20-50 years, the increase in nutrient discharged by the rivers has led to eutrophication of the coastal zones. Spring algal blooms occur frequently in Dutch coastal waters. The blooms, defined as chlorophyll a concentration ≥ 30 µg/l, are usually non-toxic, but they may be harmful because the dead algae can lead to anoxia and result in massive mussel mortality.Algal blooming is a multidisciplinary and complex problem wherehydrodynamic, chemical and biological processes take place simultaneously. The 模拟场景控制性因素光密度缺乏，不足间接的有利的协调的萌芽，发展定植模式失真的模拟有害藻类繁殖水生植物随机的决定性的详细的结构统计行为模糊逻辑技术引入制作规则表示为说明荷兰的流出物营养物泻出导致超营养作用海岸带春藻繁殖定义为叶绿素浓度无毒缺氧症大量贻贝死亡多学科的同时地Cellular Automata in Ecological and Ecohydraulics Modelling 509blooms usually occur very locally and show strong patchy dynamics. Some of the processes such as the hydrodynamics can be simulated numerically in detail, while a lot of biological mechanisms remain unclear. Besides, water quality and biological data are usually sparse and uncertain for detailed analysis.Therefore, in this research an integrated numerical and fuzzy cellular automata model was developed to predict phytoplankton biomass and hence algal blooming in the Dutch coastal waters. The numerical Delft3D-WAQ (water quality) module simulates the flow and transport conditions. The fuzzy logic module was transferredfrom the one that was developed on the basis of Noordwijk 10 data (Chen and Mynett,2004a) and was used to predict algal biomass on the basis of the computed abiotic factors. In order to take into account the spatial heterogeneity and local behaviour and to capture patchiness dynamics, a cellular automata paradigm was implemented in the developed model where the local evolution rules f are defined by fuzzy logic.The study focuses on the near shore area of the Dutch coast (Fig 7). The water depth is between 0 and 30 m, and water temperature varies from 5 to 22 ºC, and theirradiance is between 132~1700 Whm -2day -1. The concentrations of inorganic nitrogen and phosphorus are between 0.007~1.246 mg/l and 0~0.073 mg/l respectively. The biomass concentration (in chlorophyll a ) is from 0.1 ~ 90.2 µg/l.The discharge from the River Rhine at the Maassluis station (including tidal effects) is between -2744~4649 m 3/s, with a mean of 1382 m 3/s. The water is usually well mixedexcept for temporary weak stratification caused by salinity.Fig. 7. Study area of Dutch coast of the North Sea and computational grid (right)4.2 Model DevelopmentA curvilinear grid is used for the computations with Delft3-WAQ (Fig 7), totalling 1157 computational cells in the studied area. Following the computational approach as outlined in Chen (2004), the nitrogen, phosphorus and the transport are calculated first. Since the studied area is usually well-mixed with only temporary and weak stratifications, the mean water column irradiance is used.局部地局部动态数字上地稀薄的综合数值失真的浮游植物生物总量并未见有逻辑关系！模数，组件转自非生物计算空间异质性捕捉斑块动态完成近岸海域辐射度无机氮磷各自地要不是由于暂时的弱分层盐性曲线网格计数计算方法描述平均水柱辐射度510 A. Mynett and Q. ChenThe fuzzy logic model was introduced to predict algal biomass on the basis of the calculated nutrient concentrations. The membership functions of Chla and dissolvedinorganic nitrogen are shown in Fig 8. The other variables include Julian date, water temperature and concentration of dissolved inorganic phosphorus. There are in total 72 inference rules in the rule base that come both from ecologists’ experience and from the data learning process. The simulation time step for nutrients is 1 hour and that for algal biomass is 7 days, therefore a time aggregation is made before initiating the fuzzy cellular automata module. The cellular automata module is directly implemented on the topological connectivity of the curvilinear grid. The Moore neighbourhood configuration is applied in the CA model and the local evolution rules constructed by fuzzy logic techniques follow the general formula:1*1,,,(,,)t t t t i j i ji jNSf S SS ++=∑(4)where 1,t i j S +is the state of cell (i , j ) at time step t +1, *1,t i j S +is the state of the cell (i , j ) at time step t +1 which is preliminarily predicted by the fuzzy logic model, andtNS ∑are the states of the eight neighbouring cells at time step t , while f represents the local fuzzy evolution rules.Fig. 8. Membership functions of model variables and output (left: TIN; right: Chla)4.3 Results and DiscussionThe model results for dissolved inorganic nitrogen, inorganic phosphorus and chlorophyll a concentrations in 1995 are given at the left of Fig 9 while the figure at the right shows a satellite image at the peak-bloom period. The modelled peak value of chlorophyll a at station Noorwijk 10 is 48 µg/l which appeared on 30thApril. By examining the observations of 1995, it is found that the algal bloom initiated by the end of April and the first peak at station Noorwijk 10appeared on May 3rd with a chlorophyll a concentration of 58.2 µg/l. The modelled bloom timing and intensity from the fuzzy cellular automata module are seen to be somewhat close to the observations. In the spatial domain, the algal biomass is higher in the estuaries than in the coastal waters since it affected by the residual flow from the River Rhine. For comparison, the satellite image of algal bloom at the beginning of May 2003 along Dutch coast is shown in Fig 9 (right). It is difficult to quantitatively compare the modelled algae concentrations with the satellite image隶属函数溶解无机氮变量溶解无机磷时间聚合拓扑连通性配置通式预先地卫星成像峰值开始时间和强度空间范围河口残流定量对比Cellular Automata in Ecological and Ecohydraulics Modelling 511Fig. 9. Concentrations of modelled algal biomass (left) and remote sensing image (right) atpeak-bloom period.observations. However, the modelled spatial patterns are shown to be comparable with the observation. It is particularly important to point out that the application of the cellular automata paradigm has enhanced the capability of capturing patchy dynamics which can be observed from Fig 9. The case study demonstrated that this way of integrated modelling seems a promising approach indeed.5 DiscussionsIt becomes apparent from the research that in ecological modelling, cellular automata exhibit good flexibility in realising local interactions and spatial heterogeneity. In order to model practical aquatic ecosystems, it is imperative to extend the purely geometrically-determined rules to include external factors. The selection of cell size and time step is determined by the governing physical and biological processes. Based on these ideas, the CA model successfully simulated the competitive growth and succession of two underwater macrophytes in a eutrophic lake, where previous efforts using PDE models failed and had to be replaced by an individual based model.However, in aquatic ecosystems, there are usually many mechanisms remainingunclear and their statistical behaviours are difficult, if not impossible, to be defined due to limited data. Neither deterministic nor stochastic description alone is enough under such conditions. The study presented in this paper demonstrated that, although still at the initial stage, the potential to integrate different paradigms (CA and PDEs)and different techniques (fuzzy logic and numerical simulation) in ecosystem modelling. Their strength lies in the fact that (1) some processes can be simulated in detail numerically; (2) some irregular and sparse data, and the empirical knowledge from experts can be encapsulated by fuzzy logic; (3) the spatial heterogeneity, localinteractions and emergence of patchiness are captured through cellular automata.遥感空间图形证明集成建模显然的适应性实现本地互动和空间异质性实际的水生态系统势在必行的纯粹的几何确定关系营养正常的水上的生态系统潜在性不规则的，稀疏的经验主义的封装出现得到512 A. Mynett and Q. ChenReferences1.Chen, Q., Mynett, A.E., Minns, A.W., 2002. Application of cellular automata to modellingcompetitive growth of two underwater species C. aspera and P. pectinatus in Lake Veluwe. Ecological Modelling, 147: 253-2652.Chen, Q., Mynett, A.E., 2003. Effects of cell size and configuration in cellular automatabased prey-predator modelling. Simulation Modelling Practice and Theory, 11: 609-625 3.Chen, Q., 2004. Cellular Automata and Artificial Intelligence in Ecohydraulics Modelling,PhD thesis, Taylor & Francis Group Plc, ISBN: 90 5809 696 34.Marcel, S., 1999. Charophyte colonisation in shallow lakes. PhD thesis, University ofAmsterdam。