Modelling of Ecosystems by Tools from Computer Science

1. Modelling of Ecosystems by Tools from Computer Science Summer School at Czech University of Life Sciences, Prague, 16-20 September, 2013 Winfried Kurth University of Göttingen, Department Ecoinformatics, Biometrics and Forest Growth Modelling point patterns, competition and plant-herbivore interaction

2. Modelling point patterns of plant locations /* specification of an irregular stand structure in the 2D plane with random tree coordinates, but close neighbourhoods of larger individuals excluded */ module Stand; module Seedling(super.length) extends F(length, 5, 4); module Tree(super.length) extends F(length, 7, 2);

3. protected void init() [ Axiom ==> Stand; ]

4. protected void init() [ Axiom ==> Stand; ] public void make() [ Stand ==> for ((1 : n)) ( [ Translate(random(0, x_extens), random(0, y_extens), 0) Seedling(random(10, 20)) ] ); ] distribute the seedlings

5. protected void init() [ Axiom ==> Stand; ] public void make() [ Stand ==> for ((1 : n)) ( [ Translate(random(0, x_extens), random(0, y_extens), 0) Seedling(random(10, 20)) ] ); s:Seedling(h) ==> if (notOutcompeted(s) && s[Location.X] >= 0 && s[Location.Y] >= 0 && s[Location.X] <= x_extens && s[Location.Y] <= y_extens) ( Tree(factor * h + normal(0, 15)) ) else (); ] let the seedlings develop to trees if ...

6. competition rule boolean notOutcompeted(Seedling s) { return empty( (* t:Seedling, ( (t != s) && (distance(s, t) <= inhib_r) && (t[length] >= s[length]) ) *) ); }

7. parameters int n = 150; /* initial number of seedlings */ double x_extens = 500; /* extension of stand */ double y_extens = 350; double inhib_r = 35; /* distance which inhibits growth */ double factor = 4; /* proportionality of seedling and tree height */

10. Extension: consider local clustering of seedlings module Cluster(super.diameter) extends F(1, diameter, 14); module Seedling(super.length) extends F(length, 5, 4); module Tree(super.length, super.diameter) extends F(length, diameter, 2); int sd_per_cl = 20; /* number of seedlings per cluster */ int n = 8; /* number of clusters */

11. protected void init() [ Axiom ==> [ Translate(x_extens/2, y_extens/2, -1) Box(1, x_extens, y_extens).(setColor(-1)) ] /* soil */ Stand; ] visible soil

12. public void make() [ Stand ==> for ((1 : n)) ( [ Translate(random(0, x_extens), random(0, y_extens), 0) Cluster(random(60, 120)) ] ); Cluster(d) ==> for ((1 : sd_per_cl)) ( [ RH(random(0, 360)) RU(90) M(random(0, 0.5) * d) RU(-90) Seedling(random(10, 30)) ] ); s:Seedling(h) ==> if (notOutcompeted(s) && s[Location.X] >= 0 && s[Location.Y] >= 0 && s[Location.X] <= x_extens && s[Location.Y] <= y_extens) ( {double hnew = factor * h + normal(0, 20);} Tree(hnew, 0.4 * h) ) else (); ]

14. a (slightly) improved tree architecture Tree(length, diameter) ==> F(length, diameter, 2) Cone(length, length/3);

15. Tree(length, diameter) ==> F(length, diameter, 2) Cone(length, length/3);

16. Tree(length, diameter) ==> F(length, diameter, 2) Cone(length, length/3); better: instantiation rule module Tree(float length, float diameter) ==> F(length, diameter, 2) Cone(length, length/3);

17. more complex tree architecture can also be used

18. Modelling seed dispersal and plant-herbivore interaction Behaviour of plants

19. public void growPlants() [ Plant(t, r), (t > pmaxage) ==> ; Plant(t, r), (r < 0) ==> ; p:Plant, (* q:Plant *), (distance(p, q) < q[radius] && p[radius] <= q[radius]) ==> ; Plant(t, r), ((t == pgenage1 || t == pgenage2) && r >= pminrad) ==> for ((1 : (int) (pgenfac*r))) ( [ RH(random(0, 360)) RU(90) M(random(distmin, distmax)) RU(-90) Plant(0, seed_rad) ] ) Plant(t+1, r); Plant(t, r) ==> Plant(t+1, r+pgrow); ]

20. test the examples sm09_e23.rgg model for spreading (1 species) sm09_e24.rgg model for spreading (2 species) Competition is not taken into account in these examples. It is also demonstrated how population strength values can be plotted in charts during the simulation runtime.

22. Behaviour of animals

23. public void growAnimals() [ Animal(t, e), (t < 0) ==> Animal(t+1, e); /* start lag */ Animal(t, e), (e <= 0) ==> ; Animal(t, e), (e > f_e*thr) ==> [ RH(random(0, 360)) RU(90) M(shortstep) RU(-90) Animal(0, e/2 - f_e*respi) ] RH(random(0, 360)) RU(90) M(shortstep) RU(-90) Animal(0, e/2 - f_e*respi); a:Animal(t, e), (* p:Plant(u, r) *), (distance(a, p) < p[radius]) ==> RH(random(0, 360)) RU(90) M(shortstep) RU(-90) Animal(t+1, e + f_e*eat - f_e*respi) { p[radius] :-= eat; }; Animal(t, e) ==> RH(random(0, 360)) RU(90) M(longstep) RU(-90) Animal(t+1, e - f_e*respi); ]

24. protected void init() [ Axiom ==> Plant(0, seed_rad) [ RH(random(0, 360)) RU(90) M(10) RU(-90) Animal(-lag, f_e*init_e) ]; ] public void make() { growAnimals(); derive(); growPlants(); }

26. proposed model extensions: • preferred direction of seed spreading (wind) • sighted movement of animals (towards nearest plant) • more realistic reproduction (and lifespan) rules for the animals • 3-d plant architecture and animal movement • introduction of predator animals • mutation enhancing plant resistance (at cost of growth)

27. Thank you for your attention! (by Stephan Rogge)

28. Supplement Modelling above- and below-ground competition together

29. Problem: observation from ecology (A. Polle, oral communication): under optimal conditions (plenty of resources): mono-species stand gives higher yield than mixed stand under stress (lack of resources): mixed stand outperforms mono-species stand mono-species biomass multi-species stress How to model this in a simple way in GroIMP ?

30. assume 2 species, competing at two levels (e.g., above-/below-ground) • at 1st level: allelopathic interaction • i.e. species A inhibits growth of species B in its neighbourhood, • competition between A and B is stronger than between A and A or • between B and B • at 2nd level: complementarity between species A and B, • e.g. both compete for the same soil resource but at different depths •  intraspecific competition is stronger than that between A and B • Now if we decrease the availability of the soil resource (i.e., induce stress), the competition at 2nd level will gradually dominate that at 1st level; hence the advantage of the monospecies composition will vanish.

31. Objects of the model: module Stand; module Seedling(int species); module Tree(float size, int age, int species) /* size = radius */ { double nutr = satNut; } ==> if (species == 1) ( c:Cone.(setLength(3*size), setRadius(size)) {{c.setColor(0x007700);}} ) else ( Translate(0, 0, 2*size) Scale(1.0, 1.0, 2.0) s:Sphere.(setRadius(size)) {{s.setColor(0x22ee00);}} ); module Soil; module SoilPatch(int k) extends Box(3.0, 3.0, 3.0).(setColor(0xffaa88)) { float nutrAvail = startNutr; float nutrDemand, nutrTransf; int nbPlants; }

32. Competition at level 1 (allelopathy): boolean notOutcompeted(Tree a) { return empty( (* t:Tree, ((t!=a) && (( t[species] == a[species] && distance(a, t) <= inhib1 * t[size] ) || ( t[species] != a[species] && distance(a, t) <= inhib2 * t[size] )) && (t[size] >= a[size]) ) *) ); } t a const double inhib1 = 1.; /* inhibition zone: competitors of same species */ const double inhib2 = 2.; /* inhibition zone: competitors of other species */

33. Competition at level 2 (complementarity): each plant takes nutrients from all soil patches within its depletion radius species 1 from layer 1, species 2 from layer 2

34. Renewal of nutrients in a soil patch (from influx, mineralization etc.): protected void updateNutr() [ p:SoilPatch ::> { p[nutrAvail] += renewRate; p[nutrAvail] = Math.min(p[nutrAvail], patchCapacity); p[nbPlants] = 0; p[nutrTransf] = 0.; p[nutrDemand] = 0.; colorize(p); } ]

35. Checking the total demand exerted on a single soil patch: protected void checkDemand() [ t:Tree, p:SoilPatch, (inDepletionZone(p, t)) ::> { p[nutrDemand] += maxConsumption; p[nbPlants]++; } ] number of plants taking nutrients from this soil patch p

36. How to check if p is in the depletion zone of t : boolean inDepletionZone(SoilPatch p, Tree a) { return (p[k] == a[species] && projDistance(p, a) <= depletionRadius); } double projDistance(SoilPatch p, Tree a) /* distance in xy plane */ { double dx = p[Location.X] - a[Location.X]; double dy = p[Location.Y] - a[Location.Y]; return Math.sqrt(dx*dx + dy*dy); }

37. The amount of nutrient which one plant can really obtain from this soil patch is calculated and written to the soil patch attribute nutrTransf : protected void writeDemand() [ p:SoilPatch ::> { if (p[nbPlants] > 0) p[nutrTransf] = Math.min(p[nutrAvail], p[nutrDemand]) / (double) p[nbPlants]; } ]

38. Nutrient transfer from soil to plants takes place: protected void transferNutr() [ t:Tree, p:SoilPatch, (inDepletionZone(p, t)) ::> { t[nutr] += p[nutrTransf]; p[nutrAvail] -= p[nutrTransf]; } ]

39. Growth of a plant depends on its size, age, and accumulated nutrients during the last time step. All nutrients are consumed during the growth step: protected void grow() [ a:Tree(size, age, k) ==> if (notOutcompeted(a) && a[nutr] >= minNeed) ( t:Tree(growthf(size, age, a[nutr]), age+1, k) { t[nutr] = 0.; } ); ] public float growthf(double x, int age, double nutr) { double incr; if (nutr < minNeed) incr = 0.; else { if (nutr < satNut) incr = ((nutr-minNeed)/(satNut-minNeed))*maxIncr; else incr = maxIncr; } if (age < max_age) return (float) (x+incr); else return (float) x; } incr maxIncr nutr minNeed satNut

40. Initialization of the whole scene: protected void init() { step = 1; masstable.clear().setColumnKey(0, "biomass"); chart(masstable, XY_PLOT); int nbPx = (int) Math.floor((x_extens + depletionRadius) / patchDist); int nbPy = (int) Math.floor((y_extens + depletionRadius) / patchDist); for ( apply(3) ) /* 3 steps of rule application are carried out: */ [ Axiom ==> [ Translate(x_extens/2, y_extens/2, -1) /* soil surface: */ Box(0.2, x_extens, y_extens).(setColor(0xffffcc)) ] StandSoil; Stand ==> for ((1:n1)) ( [ Translate(random(0, x_extens), random(0, y_extens), 0) Seedling(1) ] ) for ((1:n2)) ( [ Translate(random(0, x_extens), random(0, y_extens), 0) Seedling(2) ] ); Soil ==> for (int i: (1:nbPx)) ( for (int j: (1:nbPy)) ( for (int k: (1:2)) /* two layers of patches */ ([ Translate(-0.5*depletionRadius + i*patchDist, -0.5*depletionRadius + j*patchDist, -k*layerDist) p:SoilPatch(k) { colorize(p); } ]) ) ); s:Seedling(k) ==> if (s[Location.X] >= 0 && s[Location.Y] >= 0 && s[Location.X] <= x_extens && s[Location.Y] <= y_extens) ( Tree(random(1.0, 5.0), 1, k) ); /* size = random, age = 1, species = k */ ] }

41. Colorizing the patches for visualization of nutrient content: void colorize(SoilPatch p) { double a = p[nutrAvail]; if (a >= 0.875*patchCapacity) p.setColor(0xff4444); else if (a >= 0.75*patchCapacity) p.setColor(0xff7744); else if (a >= 0.625*patchCapacity) p.setColor(0xffaa44); else if (a >= 0.5*patchCapacity) p.setColor(0xffdd44); else if (a >= 0.375*patchCapacity) p.setColor(0xffff55); else if (a >= 0.25*patchCapacity) p.setColor(0xffff88); else if (a >= 0.125*patchCapacity) p.setColor(0xffffbb); else p.setColor(0xffffff); } (red  yellow  white)

42. Main rule block which controls the order of rule applications and adds a biomass entry after each growth step for a chart: public void run() { updateNutr(); checkDemand(); writeDemand(); transferNutr(); grow(); masstable.addRow().set(0, sum( biomass( (* t:Tree *) ))); println(step); step++; } public float biomass(Tree a) /* assumed as proportional to basal area */ { return (float) (Math.PI * a[size] * a[size]); }

43. for completeness: parameters of the model and variable declarations const int n1 = 200; /* nb. of trees of species 1 */ const int n2 = 200; /* nb. of trees of species 2 */ const int max_age = 40; /* number of growth steps */ const double x_extens = 500.; const double y_extens = 300.; const double inhib1 = 1.; /* inhibition zone: competitors of same species */ const double inhib2 = 2.; /* inhibition zone: competitors of other species */ const double renewRate = 20.0; /* let this vary! */ const double patchCapacity = 50.0; const double startNutr = 50.0; const double maxConsumption = 0.4; /* max. consumption of a tree from a single soil patch; 0.4 */ const double patchDist = 15.; const double depletionRadius = 40.; const double layerDist = 5.; const double minNeed = 4.; /* necessary min. amount of nutrients for survival */ const double satNut = 10.; /* saturation point of nutrient response function */ const double maxIncr = 0.7; /* max. tree growth rate */ int step; /* step counter */ const DatasetRef masstable = new DatasetRef("stand biomass"); /* table for chart */ example file div_compet2b.rgg

45. Resulting stand biomass under different parameterizations (unstressed = plenty of nutrients by high renewal rate in the soil; stressed = competition for nutrients is severely limiting)