Ecological risk assessment of genetically modified crops based on cellular automata modeling

Biotechnology Advances 27 (2009) 1132–1136

Contents lists available at ScienceDirect

Biotechnology Advances j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / b i o t e c h a d v

Ecological risk assessment of genetically modified crops based on cellular automata modeling Jun Yang a,b, Zhi-Rui Wang c, De-Li Yang b, Qing Yang a, Jun Yan a, Ming-Feng He c,⁎ a b c

Department of Bioscience and Biotechnology, Dalian University of Technology, Linggong Road No.2, Dalian, 116024, PR China School of Management, Dalian University of Technology, Linggong Road No.2, Dalian, 116024, PR China Department of Applied Mathematics, Dalian University of Technology, Linggong Road No.2, Dalian,116024, PR China

a r t i c l e

i n f o

Available online 27 May 2009 Keywords: Genetically modified (GM) crops Ecological impact Risk assessment Cellular automata Simulation

a b s t r a c t The assessment of ecological risk in genetically modified (GM) biological systems is critically important for decision-making and public acceptance. Cellular automata (CA) provide a potential modeling and simulation framework for representing relationships and interspecies interactions both temporally and spatially. In this paper, a simple subsystem contains only four species: crop, target pest, non-target pest and enemy insect, and a three layer arrangement of L × L stochastic cellular automata with a periodic boundary were established. The simulation of this simplified system showed abundant and sufficient complexity in population assembly and densities, suggesting a prospective application in ecological risk assessment of GM crops. © 2009 Elsevier Inc. All rights reserved.

1. Introduction Despite increasing adoption (James, 2008) and numerous beneficial promises (Huang et al., 2002; Chapotin & Wolt, 2007), there is a great deal of concern about the potential undesirable impacts of genetically modified (GM) crops on the environment, health and economics (Kuiper et al., 2001; Conner et al., 2003; Belcher et al., 2005; Zelaya et al., 2007; Schubert, 2008; Finamore et al., 2008). Risk assessment provides a critically important control procedure for decision-making and to promote the public acceptance of GM crops. However, the prediction of risk associated with genetically modified biological systems is not easy due to the complexity of natural environments and ecosystems, and the need for risk to be considered at multiple spatial or temporal scales (Hill, 2005). Ecological risk assessment (ERA) provides a robust process for the management of environmental risks (Dale et al., 2008). Although there are different regulatory frameworks worldwide (Sparrow, 2009; Paoletti et al., 2008), ERA is generally characterized by problem formulation, analysis of exposure and exposure–response relationships, risk characterization and communication (EFSA, 2006, 2008). The most common approaches used are unstructured brainstorming, checklists, Hierarchical Holographic Modeling (Hayes et al., 2004), tier tests (Garcia-Alonso et al., 2006), and integrated assessment frameworks (Cockburn, 2002; Varzakas et al., 2007; Vergragt and Brown,

⁎ Corresponding author. E-mail address: [email protected] (M.-F. He). 0734-9750/$ – see front matter © 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.biotechadv.2009.05.024

2008). However, ERA has been challenged regarding quantitative uncertainty analysis of the dynamic ecosystem (Hayes et al., 2004), large scale temporal and spatial analysis (Hill, 2005; Belcher et al., 2005; Munro, 2008), obstacles of standardization deficiency (Hill, 2005) and data insufficiency (Raybould, 2006; Lang et al., 2007). The development of mathematical models is of interest in ecological assessment, since models allow us to understand and predict phenomena at a variety of spatial and temporal scales. Cellular Automata (CA) is a discrete, spatially explicit extended dynamic system that can react to external stimuli (input) with a state transition and an output value according to given functions, at discrete time intervals (Wolfram, 1984; He et al., 2002). Because CA models are not equation based, this allows for the direct consideration of knowledge that is not necessarily restricted to a large amount of data (Jeltsch et al., 1996) and is particularly useful in consideration of complex systems (Wolfram, 1984). CA models have shown promising application and have been used to simulate spatial dynamics of GM crops, including gene flow (Linacre & Ades, 2004), the evolution of resistance to Bacillus thuringiensis (Bt) toxin in transgenic maize (Peck et al., 1999) and cotton (Tyutyunov et al., 2007), the spatial patterns of GM crop contamination (Belcher et al., 2005; Munro, 2008), and the spatial configuration of refuges for resistance management of Bt transgenic crops (Vacher et al., 2003; Cerda & Wright, 2004). The temporal and spatial aspect of the CA model provides advantages for simulating population distribution and interactions, which makes the CA model well suited for the purpose of ERA of GM crops. GM Bt crops are widely planted (James, 2008). A central concern in regulating Bt crops is the risk of insects evolving resistance to the Bt toxins. However, recent investigations have shown that secondary

J. Yang et al. / Biotechnology Advances 27 (2009) 1132–1136

pests, such as mirids, which are not targeted by the Bt toxin and used to be controlled by broad-spectrum pesticides have already become a new threat in Bt crops (Wang et al., 2006; Wang et al., 2008; Wu et al., 2008; Naranjo et al., 2008). Because of the high mobility, broad host range, critical damage and poorly understood population dynamics of mirid pests, their monitoring and management are often problematic (Lu et al., 2008). Mathematical models may provide a possible solution in predicting the population dynamics of these plant bugs. In this research, we focus on the modeling of population dynamics in a simple subsystem containing plants (cotton), target pests (cotton bollworm), non-target secondary pests (mirid) and enemy insects by Cellular Automata. We expect the application of this model will help us understand the potential evolution of this system at certain spatial and temporal scales, and provide information that is helpful in ecological risk assessment and pest management for transgenic Bt crops.

2. The model 2.1. Model description We define a three-layer L × L stochastic cellular automata (Fig. 1) with a boundary period where L is an even number. For any lattice, denoted as a matrix M, in which M(m,n) represents the location of each cell by line m and row n. Each gird type was defined as γ, γ∈ (0,1). The defined species D and blanks were randomly placed on the lattice layer M1; species B, C and blanks were randomly placed on layer M2. Each grid is occupied, at most, by an individual. Species A is the only species which has been quantified and was placed at each lattice point on layer M3. A contiguous pair of grid α and β change with the probability of P, either (α,β) to (α',β'), or remain unchanged, depends on whether the gird is occupied by an individual or not. The relationships between these four species were described and evolved through long term scales.

Fig. 1. Scheme of a three-layer of L × L stochastic cellular automata.

1133

Fig. 2. The relation between four species.

2.2. The definition of individuals and simulation steps In order to describe the relationships between the four species, we define A, B, C, D as plant, target pest, non-target secondary pest, and enemy insect respectively, and denote “→” as predation. We then have the following predation relations (Fig. 2): B →A, C →A, D →B, D →C. As the resource supplier, only species A was quantified in this model, and individual variation was ignored. The quantity of A at each grid was denoted by nA(i,j) ≥ 0. Each individual of species B, C, and D was given detailed characteristics, which were described by the Penna asexual model (Penna, 1995). We define hBi ,hCi ,hD i for the hunger level of individuals of species B C D , hmax , hmax for their maximum endurable B, C, and D; and hmax hunger. If the hunger level of an individual is larger than its maximum endurable hunger, it will die of starvation. We define the mature age of each species B, C, and D as r B,rC,r D. If an individual reaches its mature age, it will generate an offspring. We define the probability of predation occurring as PBA,PCA,PDB,PDC. If predation succeeds, the predator's hunger will be reduced by hBA, hCA, hDB, hDC, and the prey will die and vacate the occupied grid. If the predation fails, the prey will escape and hunger of the predating individual will increase by UB, UC, UD. As described above (see 2.1, Fig. 1), species A is quantified, as defined by nA(i,j).We provide definitions dB and dC for the reduction of n A(i,j), if it was predated by species B or C. When t = 0, the quantity of species A is maximum, denoted by n A(i,j) = Cmax; the individuals of species B, C and D and blanks in layers M2 and M3were set randomly with a probability of P = 0.5. When given a random time, t, we can then follow the three-step running of this model as: First, we simulate stochastic cellular automata of each individual horizontally at layer M1, M2. For each contiguous pair of grid a and b, only when the following two conditions satisfied will the states changed with a probability of P= 1: if an individual has not reached its mature age, the states changed from (1,0) to (0,1), and (0,1) to (1,0); if an individual has reached its mature age, the states changed from (1,0) to (1,1), and (0,1) to (1,1),where 1 represents the new born. Then, we stimulate stochastic cellular automata of each individual vertically, and denote M1(i,j,t), M2(i,j,t), M3(i,j,t). If there is an individual on M1(i,j,t), it is labeled as D, otherwise as 0. Any individual on M2(i,j,t) is labeled as B and C, if present, or 0. The quantity of A on M3(i,j,t) is denoted as nA(i,j,t). When D predation on B succeeds, the starvation of D is reduced by hDB; when D predation on C succeeds, the starvation of D is reduced by hDC; when B predation on A succeeds, the starvation of B is reduced by hBA; when C predation on A succeeds, the starvation of C is reduced by hCA. Finally, we check the individuals' survival state according to Penna model. The increase in starvation is measured by UA, UB, UC, UD. If individual dies, the grid becomes vacant. Then, the system will run the next simulation cycle.

1134

J. Yang et al. / Biotechnology Advances 27 (2009) 1132–1136

2.3. Simulation results The simulation was first run with the main parameters as L = 500, t = 2000, P DB = P DC > 0.5. We focus on the repeated simulation of interspecies' interactions based on the full probability of B → A and C → A (PBA, P CA). The representative results from the above simulation are shown in Fig. 3. Although we define a higher but equal feeding probability for species D to B and to C, the population structure in this subsystem varied in a very dynamic way after 2000 simulation cycles. The population size of species A was influenced significantly by preferential feeding by pests B and C (Fig. 3a). Feeding preference is also a key factor in the population dynamics of species B and C, which indicates a competitive relationship between the plant pests. They actually shared very modest population spaces (Fig. 3b,c). This result also indicates that the reduced feeding preference of the target pest (equal to the repression of the target pest by the Bt toxin) will definitely result in the increase in population of the non-target pest. Another simulation was conducted with the main parameters at L = 500, t = 2000, and PBA > PCA ≥ 0.5. We focused on the repeated simulation of interspecies' interactions based on the full probability of D → B and D → C (PDB, PDC). The results from this simulation are shown in Fig. 4. In this simulation, the feeding preference of D (on this occasion, D also represents other selective pressures, such as insecticide sprays and Bt toxins.) on pests B and C showed its power in controlling the community assembly in this simple system (Fig. 4a,b,c). The feeding capacity of the enemy insects dramatically

affected the species composition and the population size of the plant pests (Fig. 4b, c), contributing to a higher yield of the plants (Fig. 4a). This simulation also indicates that the repression or removal of a primary target pest will gradually allow for a possible outbreak of non-target secondary pests. 3. Concluding remarks The assessment of ecological risk in genetically modified biological systems is critically important for decision-making and public acceptance. As a discrete, spatially explicit extended model, CA showed a high potential for application in representing the relationships and interactions within a simple system, and generated varied dynamics of population assembly. Although the conceptual plant-pests-enemy CA model we created in this research is simplified showing only the basic components of a natural subsystem, the simulation results showed abundant and sufficient complexity in population assembly and population densities. Further development of this model with practical data acquisition and input may provide more informative details on the ecological consequences of transgenic crops. Acknowledgments The authors acknowledge the financial support provided by the National R&D Project of Transgenic Crops of Ministry of Science and Technology of China (JY03-B-18-02), Municipal Science and Technol-

Fig. 3. Simulation of population dynamics at full probability of B → A and C → A (PBA, PCA). (a) The average population percentage of species A. (b), (c) and (d) The population dynamics of species B, C and D.

J. Yang et al. / Biotechnology Advances 27 (2009) 1132–1136

1135

Fig. 4. Simulation of population dynamics at full probability of D → B and D → C (PDB,PDC). (a) The average percentage of A. (b), (c) and (d) The population dynamics of B, C, D.

ogy Project of Dalian (2007B10NC137), Key Laboratory of Industrial Ecology and Environmental Engineering (Dalian, China), and Key Laboratory of Natural Pesticide & Chemical Biology (Wuhan, China). References Belcher K, Nolan J, Phillips PWB. Genetically modified crops and agricultural landscapes: spatial patterns of contamination. Ecol Econ 2005;53:387–401. Cerda H, Wright DJ. Modeling the spatial and temporal location of refugia to manage resistance in Bt transgenic crops. Agri Eco Environ 2004;102:163–74. Chapotin SM, Wolt JD. Genetically modified crops for the bioeconomy: meeting public and regulatory expectations. Transgenic Res 2007;16:675–88. Conner AJ, Glare TR, Nap JP. The release of genetically modified crops into the environment Part II. Overview of ecological risk assessment. Plant J 2003;33:19–46. Cockburn A. Assuring the safety of genetically modified (GM) foods: the importance of a holistic, integrative approach. J Biotechnol 2002;98:79-106. Dale VH, Biddinger GR, Newman MC, Oris JT, Suter GW, Thompson T, et al. Enhancing the ecological risk assessment process. Integr Environ Assess Manag 2008;4:306–13. European Food Safety Authority (EFSA). Guidance document of the scientific panel on genetically modified organisms for the risk assessment of genetically modified plants and derived food and feed. EFSA J 2006;99:1–100. European Food Safety Authority (EFSA). Report of the EFSA GMO Panel Working Group on animal feeding trials. Safety and nutritional assessment of GM plants and derived food and feed: the role of animal feeding trials. Food Chem Toxicol 2008;46:S2–S70. Finamore A, Roselli M, Britti S, Monastra G, Ambra R, Turrini A, et al. Intestinal and peripheral immune response to MON810 maize ingestion in weaning and old mice. J Agric food Chem 2008;56:11533–9. Garcia-Alonso M, Jacobs E, Raybould A, Nickson TE, Sowig P, Willekens H, et al. A tiered system for assessing the risk of genetically modified plants to non-target organisms. Environ Biosafety Res 2006;5:57–65.

Hayes KR, Gregg PC, Gupta VVSR, Jessop R, Lonsdale M, Sindel B, et al. Identifying hazards in complex ecological systems. Part 3: Hierarchical Holographic Model for herbicide tolerant oilseed rape. Environ Biosafety Res 2004;3:1–20. He MF, Lin J, Jiang H, Liu X. The two populations' cellular automata model with predation based on the Penna model. Physica A 2002;312:243–50. Hill RA. Conceptualizing risk assessment methodology for genetically modified organisms. Environ Biosafety Res 2005;4:67–70. Huang J, Pray C, Rozelle S. Enhancing the crops to feed the poor. Nature 2002;418:678–84. James C. Global Status of Commercialized Biotech/GM Crops: 2008. ISAAA Brief No. 39. Ithaca, NY: ISAAA; 2008. Jeltsch F, Milton SJ, Dean WRJ, VanRooyen N. Tree spacing and coexistence in semiarid savannas. J Ecol 1996;84:583–95. Kuiper HA, Kleter GA, Noteborn HPJM, Kok EJ. Assessment of the food safety issues related to genetically modified foods. Plant J 2001;27:503–28. Lang A, Lauber E, Darvas B. Early-tier tests insufficient for GMO risk assessment. Nat Biotechnol 2007;25:35–6. Linacre NA, Ades PK. Estimating isolation distances for genetically modified trees in plantation forestry. Ecol Model 2004;179:247–57. Lu YH, Qiu F, Feng HQ, Li HB, Yang ZC, Wyckhuys KAG, et al. Species composition and seasonal abundance of pestiferous plant bugs (Hemiptera: Miridae) on Bt cotton in China. Crop Prot 2008;27:465–72. Munro A. The spatial impact of genetically modified crops. Ecol Econ 2008;67:658–66. Naranjo SE, Ruberson JR, Sharma HC, Wilson L, Wu KM. The present and future role of insect-resistant GM crops in cotton IPM. In: Romeis J, Shelton AM, Kennedy GG, editors. ntegration of insect-resistant GM crops within IPM programs. Dordrecht, the Netherlands: Springer; 2008. p. 159–94. Paoletti C, Flamm E, Yan W, Meek S, Renckens S, Fellous M, et al. GMO risk assessment around the world: some examples. Trends Food Sci Technol 2008;19:S70–8. Peck SL, Gould F, Ellner SP. Spread of resistance in spatially extended regions of transgenic cotton: Implications for management of Heliothis virescens (Lepidoptera: Noctuidae). J Econ Entomol 1999;92:1–16. Penna TJP. A bit-string model for biological aging. J Stat Phys 1995;78:1629–33. Raybould A. Problem formulation and hypothesis testing for environmental risk assessments of genetically modified crops. Environ Biosafety Res 2006;5:119–25.

1136

J. Yang et al. / Biotechnology Advances 27 (2009) 1132–1136

Schubert DR. The problem with nutritionally enhanced plants. J Med Food 2008;11:601–5. Sparrow PAC. GM Risk assessment. In: Jones HD, Shewry PR. editors. Methods in Molecular Biology, Transgenic Wheat, Barley and Oats. Human Press, LLC; 2009. Tyutyunov YV, Zhadanovskaya EA, Arditi R, Medvinsky AB. A spatial model of the evolution of pest resistance to a transgenic insecticidal crop: European corn borer on Bt maize. Biophysics 2007;52:52–67. Vacher C, Bourguet D, Rousset F, Chevillon C, Hochberg ME. Modelling the spatial configuration of refuges for a sustainable control of pests: a case study of Bt cotton. J Evol Biol 2003;16:378–87. Varzakas TH, Chryssochoidis G, Argyropoulos D. Approaches in the risk assessment of genetically modified foods by the Hellenic Food Safety Authority. Food Chem Toxicol 2007;45:530–42.

Vergragt PJ, Brown HS. Genetic engineering in agriculture: new approaches for risk management through sustainability reporting. Technol Forecast Soc Change 2008;75:783–98. Wang S, Just DR, Pinstrup AP. Damage from secondary pests and the need for refuge in China. In: Just RE, Alston JM, Zilberman D, editors. Regulating Agricultural Biotechnology: Economics and Policy. New York: Springer; 2006. p. 625–37. Wang S, Just DR, Pinstrup AP. Bt-cotton and secondary pests. Int J Biotechnol 2008;10:113–21. Wu KM, Lu YH, Feng HQ, Jiang YY. Suppression of cotton bollworm in multiple crops in China in areas with Bt toxin-containing cotton. Science 2008;321:1676–8. Wolfram S. Universality and complexity in cellular automata. Physica D 1984;10:1–35. Zelaya IA, Owen MDK, VanGessel MJ. Transfer of glyphosate resistance: evidence of hybridization in Conyza (Asteraceae). Am J Botany 2007;94:660–73.