Analytical Trophodynamics Applied to Modeling Forest Dynamics with Carbon Cycling

Abstract

Models based on analytical trophodynamics (AT) method have provided an analytical framework for modeling in ecology, including the dynamical flux of nutrients present in the soil for a fixed region. Dynamics occurring concurrently in different time scales are modeled. Through a mathematical treatment of the elements of both biotic and abiotic factors, is established  stability and conservation laws for growing trajectories, whose solutions of the second-order differential systems equations known as Volterra–Hamilton systems. All solutions trajectories obtained to follow the biological principles of energy conservation. The tensors of AT were computed with the computational algebraic package FINSLER. Moreover, in this chapter, we present an overview of the last results and actual status of research in this area.

1. Introduction

Forest dynamics are closely related to the effects of climate change, biodiversity, air temperature regulation, precipitation patterns, carbon cycling, and fixation of essential nutrients for organic life such as nitrogen, phosphorus, and potassium [1]. Understanding the cycling of nutrients and carbon is essential to comprehend this process [2]. Such cycles present in forest ecosystems are influenced by biotic and abiotic effects, including ecological interactions and adaptations to environmental external changes [3, 4]. In [5], experimental results for species in Tropical Forests, show the role of the understory in the dynamics of interaction between different species of underground plants and microorganisms. In particular, it has already been shown that ectomycorrhiza fungi are responsible for the hyphae pathways that mediate transports in forest soil, increasing the rate of nutrients uptake, weathering of soil minerals, and also influencing the bacterial community in the rhizosphere [6]. Recent research [4], investigating the transfer of carbon between trees in a fixed area, under the effect of a specific microclimate, highlights the role of the rhizosphere in nutrient dynamics. Mycorrhizal networks play the role of mediators of interactions between trees and plants surrounding the understory [7, 8, 9].

The transfer of carbon and other nutrients between plants in the same place has been investigated, taking into account the mycorrhizal networks, which provide a means for the sending and receiving of nutrients to occur. Understanding quantitative aspects of the direction and magnitude, among other factors, of the transfer process directly affects the understanding of mechanisms related to the interaction and function of these plants [4, 6, 10].

Recent results [11] show that mycorrhizal fungi influence the coexistence and diversity of plants, being mediators of the symbiosis process of intra- and interspecific competition of plant communities. The different types of mycorrhizal fungi affect the level of nutrient sharing between neighboring plants and, consequently, their availability, acting as a regulator and stabilizer of the interaction mechanisms that lead to the coexistence of species.

To evaluate the concentration of nutrients present in the subsoil of regions where the different species are found, a long-term analysis under the effects of external agents, such as natural disasters and anthropic interference, is essential [12]. The long-term facilitation-type interactions have been shown to be more influential in the dynamics than the competition. This is characterized by the positive effect in which at least one of the local species benefits from interactions. A direct implication of this result is that plants with higher biomass may not prevail in the long term [13, 14].

Models based on the analytical trophodynamics (AT) method [15, 16] provide a new theoretical approach for a mathematical analysis of the carbon cycle dynamics and of other nutrients present in the soil, taking into account short-term and long-term as competition and facilitation, respectively. In addition, abiotic disruptive factors can also be integrated. AT has made it possible to model complex phenomena such as forest dynamics [5, 9, 17, 18]. Recent AT results describe the interaction between trees and the understory, including litter [19], succession forest [20], and the environmental impact of the transgenic crop [21].

Models based on TA can include external environmental effects in the local dynamics by analogy with continuum mechanics. In this approach, geometric quantities described through Finsler geometry and systems of second-order differential equations provide mathematical objects that allow the analysis of these dynamics. The algebraic manipulation of the geometric elements modeled is usually performed with the FINSLER package [22].

The effect of external forces on the local dynamics of plants affects their forms of interaction and functionality, causing stability scenarios to be altered. The eco-strain model [20] makes it possible to investigate the external effects of the surrounding space on the forest dynamics. Thus, transformations like that are important to deform the Huxley allometric space, in which the known classical ecological interactions (symbiosis, competition, and parasitism) occur so that in the new space there are optimal trajectories that conserve Medawar growth energy [16], being described by analytical expressions that modeling this mechanism.

The present work provides a review of the article [23], detailing biological and mathematical aspects of the constructed model, with the TA method. The mathematical model built using this approach took into account the concept of the potential of the mycorrhizal network [24] so that the solutions of the proposed model allowed us to describe the dynamics of carbon cycling for two species of trees in the British forest, corroborating experimental results. The trajectories of the proposed system of differential equations describe the dynamics of species in a production space disturbed by external agents. The results achieved in this model can be statistically tested at later times since it is expected that the structural behavior of the dynamics will continue to be valid for particular cases, in which the coefficients assume numerical values measured through field experiments.

In Section 1, we present an updated overview of the importance of the type of modeling. Sections 2 and 3 introduce the theoretical framework of the proposed approach. Examples of modeling applications with real data from biological experiments are shown in Section 4. We end with a description of the next steps to continue this research.

2. Volterra: Hamilton systems and Finsler and Wagner spaces

Given a coordinate chart $\left(x^i,y^i\right)$ in the base space $M$ of a slit tangent bundle ${\mathit{TM}}^{\ast }$ consider the Volterra–Hamilton [25] system defined by²

$$\begin{array}{l}\frac{{\mathit{dx}}^i}{\mathit{dt}}=k_{\left(i\right)}{y}^i,\\ \frac{{\mathit{dy}}^i}{\mathit{dt}}=G_{\mathit{jk}}^i{y}^j{y}^k+r_j^i{y}^j+e^i,i,j,k=1,\dots ,n\end{array}$$E1

where $k_{\left(i\right)}$³ is rate growth per capita, $G_{\mathit{jk}}^i\left(x,y\right)$ are positively homogeneous of degree zero in $y^i$ representing the interaction between populations $y^j$ and $y^k,r_j^i$⁴ and $e^i$ are, respectively, how the populations grow and react to environmental influence. The variable $x$ is the variable of production of Volterra [26]. It measures the quantity produced by a population $N^i=\frac{{\mathit{dx}}^i}{\mathit{dt}}$ along an interval of time.

There are three possibilities for the interaction coefficients: (i) are all constants; (ii) depends explicitly of $x^i$; or (iii) depends on the ratio $\frac{y^i}{y^j}$ [16]. When (iii) occurs, the effects of social interaction [27] are included in modeling. Furthermore, a condition assumed in this work is that the population are proportional, that is, the ratios $\frac{N^i}{N^j}$ are constants. This is known as how uniformity condition.

Given an initial condition, solutions of the second order system (SODE) (1) provide trajectories that have the property of conserve growth energy [28] along the solutions curves in the correspondent production process space. In [16] is studied a two-dimensional case of this SODE model with classical ecological interactions between two populations. In this case, the $G_{\mathit{jk}}^i$ will be constant along the time.

For the case that the coefficients $G_{\mathit{jk}}^i$ are constant these equations define the coefficients of the Berwald connection [15, 29, 30, 31]. In this case, three of the eight SODEs, are elastic transformation and are denoted by $A_{\mathit{jk}}^i$. They are the three cases are the Minkowski space $\left(M,F\right)$, where $M$ is as before, and $F$ is the Finsler norm such that with $\mathit{ds}=F\left(y\right)$, and the metric tensor⁵

$$g_{\mathit{ij}}\left(y\right)=\frac{1}{2}{\dot{\partial }}_i{\dot{\partial }}_j{A}_{\alpha }^2,$$E2

for

$$A_{\alpha }=exp\left[\varphi \left(x\right)\right]F\left(y\right),$$E3

where there are three possibilities for $F$, given by $F=F_{\alpha }\left(y\right)$, $\varphi ={\varphi }_{\alpha }\left(x\right)$, $\alpha \in \left\{\mathrm{1,2,3}\right\}$, as stated in [20].

The ideal case of Huxley allometric space is one of these, and your geodesics are given by

$$\frac{d^2{x}^i}{{\mathit{ds}}^2}=0,$$E4

for a parameter $s$ is known as a natural parameter [16]. They are the perfect case of Huxley allometric space. Consequently, the solutions are straight lines $y^i=a^{i}s+b^i$, for arbitrary real constants of integration $a^i$ and $b^i$.

As in forest ecosystems, it is common to have processes occurring on different time scales [7, 14], in the modeling developed with TA different scales are combined. While the parameter $t$ indicates the real-time scale, the parameters $s$ and $p$ – introduced ahead – denote other time scales, in which forest dynamics phenomena studied here have patterns that preserve quantities of interest, obeying the principle of conservation of the energy in the allometric space. The Finslerian geodesics given by the metric tensor (2) preserve energy growth along trajectories in production process space.

Here transformations of the perfect case are reversible in the sense that one can be obtained from the other by a non-singular change of coordinates $x^i$ into ${\overline{x}}^i$,

$$d{\overline{x}}^i=J_j^i\left(x\right){\mathit{dx}}^j,$$E5

where the Jacobian $J_j^i=\left(\frac{\partial {\overline{x}}^i}{\partial x^j}\right)$ depends only on the Volterra production variable $x^i$.

The other five, of all eight SODEs mentioned above, are plastic deformation and denoted by $A_{\mathit{jk}}^i+W_{\mathit{jk}}^i$ where

$$W_{\mathit{jk}}^i={\delta }_j^i{\partial }_k\left[\sigma \left(x\right)\right],$$E6

for a scalar function $\sigma \left(x\right)$ defined over $M$, and ${\partial }_k\left[\sigma \left(x\right)\right]={\sigma }_k\left(x\right)$ is your gradient.

According to the TA approach, after a heterochronic transformation $\left\{\sigma \left(x\right),c^i(x,y)\right\}$, we obtain the following perturbed SODE from (1):

$$\begin{array}{l}\frac{d^2{x}^i}{{\mathit{dp}}^2}+\left[A_{\mathit{jk}}^i+{\delta }_j^i{\sigma }_k\right]\frac{{\mathit{dx}}^j}{\mathit{dp}}\frac{{\mathit{dx}}^k}{\mathit{dp}}=c^i,\\ \frac{\mathit{dp}}{\mathit{ds}}=Qexp\left[-{\int }_{\gamma }\sigma \left(r\right)\mathit{dr}\right],\\ \mathit{ds}=exp\left(\mathit{\lambda t}\right)\mathit{dt},\end{array}$$E7

where $Q$ is an arbitrary constant and $\gamma$ is a solution curve. When $c^i\ne 0$, then the environmental effects act, such that the populations dynamics react through a time-sequence of growth along their trajectories [32].

Now the solutions of this new SODE (7), the autoparallels of a constant Wagner connection, are not geodesics [33]. In [34], it was demonstrated that the energy along these paths is conserved, and they provide a condition for optimality accordly division of labour principle [16, 28, 35].

Although model (7) includes the effect of environment in dynamics, it does not explain how soil carbon is periodically cycled and how it occurs in the case of interactions of paper birch and Douglas-fir trees [6]. The energy flux and mycorrhizal networks was introduced in the model by consider a proper form of the anholonomic frames $B_j^i\left(x,\mathit{dx}\right)$, depend on $x^i$ and ${\mathit{dx}}^i$, that describe the plastic deformation. Through this element, defined below, the new version of SODE (7) will provide solutions that model the nutrient cycle as expected.

External events can have a high impact on forest dynamics, directly affecting soil organic matter, and leading most microorganisms to suffer drastic reductions in their population due to high temperatures [36]. In addition to the direct effects of short-term wildfires affecting the soil microbial communities, there are indirect abiotic effects. In the long term, there are changes in both microbial and plant communities, chemical changes at different soil levels, alterations in the microclime, reduction of overstory and understory plant diversity, soil microflora, and others [37]. In turn, such effects imply patterns of forest succession such that the vegetation cover present as well as the dynamics of this ecosystem are altered [38]. Although there are approaches based on Finsler geometries, different from TA, for modeling this aspect [39, 40], understanding the long-term effects of wildfires on forest dynamics is still a research area with many challenges [41, 42, 43].

3. Randers space

For investigating the effect of surround local environment, not considered by the initial Eco-strain model [20], in [23], it was assumed that there was an increase in the total intake of nutrients per unit of time, according to the biological event [44]. Mathematically, this corresponds to assume the following Finsler metric functional

$${\mathcal{F}}_{\alpha }=exp\left[\sigma \left(x\right)\right]A_{\alpha }+b_i\left(x\right)y^i,$$E8

and its norm ${\mathit{dS}}_{\alpha }={\mathcal{F}}_{\alpha }\left(x,\mathit{dx}\right)$. The metric tensor in this case is denoted by $h_{\mathit{\alpha \beta }}\left(x,y\right)$. Using parameter $p$ of Wagner SODE, we can obtain [23, 45] your geodesics,

$$\begin{array}{l}\frac{{\mathit{dy}}^i}{\mathit{dp}}+\left[A_{\mathit{jk}}^i+{\delta }_j^i{\sigma }_k\left(x\right)+{\delta }_k^i{\sigma }_j\left(x\right)-g_{\mathit{jk}}{\sigma }^i\right]y^j{y}^k=g^{\mathit{im}}\left(y\right)\left({\delta }_m{b}_l-{\delta }_l{b}_m\right)y^l,\\ y^i=\frac{{\mathit{dx}}^i}{\mathit{dp}},{\delta }_m={\partial }_m-A_m^r{\dot{\partial }}_r.\end{array}$$E9

The ${\delta }_m$ is known as Berwald delta derivative⁶ and $g^{\mathit{im}}$ denotes the inverse of metric tensor.

Another geometric element for the model proposed in [23] is the introduction of the infectivity potential represented by the real positive function $V\left(x\right)$. Considering how a function splits off the vector field $\left(b_i\right)$ into $b_i\left(x\right)={\partial }_{i}V\left(x\right)+{\tilde{b}}_i\left(x\right)$, the potential $V\left(x\right)$ is defined so that increase the total intake per unit of time because

$${\mathit{dS}}_{\alpha }=exp\left[\sigma \left(x\right)\right]A_{\alpha }+\mathit{dV},$$E10

where $\mathit{dV}>0$.

In [23], the effect of a strained environment pressuring the forest dynamics was modeled by the Eco-strain model [20] where the Berwald frame $B_i^{\alpha }\left(x,\mathit{dx}\right)$ has been used, and additionally by the Holland’s frame $U_i^{\alpha }\left(x,\mathit{dx}\right)$.

From the original allometric space described in terms of the metric tensor (2), the relaxed state${\delta }_{\mathit{\alpha \beta }}$ and strained state$g_{\mathit{ij}}$ are obtained so that

$$\begin{array}{l}{\delta }_{\mathit{\alpha \beta }}\left(x,y\right)=g_{\mathit{ij}}\left(x,y\right)B_{\alpha }^i{B}_{\beta }^j,\\ g_{\mathit{ij}}\left(x,y\right)={\delta }_{\mathit{\alpha \beta }}{B}_i^{\alpha }{B}_j^{\beta }.\end{array}$$E11

When the frame is a Jacobian matrix it is called holonomic (or elastic) and indicates a coordinate change that is reversible. Already when this does not occur, the frame is termed anholonomic (or plastic), and it is non-reversible.

Additionally, an already deformed space can be deformed again through a composition of frames. From the Holland’s frame for Randers space, denoted by $U_i^{\alpha }$ and defined by

$$\begin{array}{l}{g}_{\mathit{ij}}\left(x,y\right)=U_i^{\alpha }{U}_j^{\beta }{h}_{\mathit{\alpha \beta }},\\ h_{\mathit{\alpha \beta }}\left(x,y\right)=g_{\mathit{ij}}{U}_{\alpha }^i{U}_{\beta }^j.\end{array}$$E12

obtain one of the plastic deformation of the Euclidean allometric space with a non-vanishing curl.

4. Applying the modeling in a concrete case

In [23], a model for carbon and nutrient cycling including understory and underground mycorrhizal networks was presented. This model extends the Eco-strain model [16]. Now, the results proposed allow to describe the periodicity of carbon cycling. The effect of mycorrhizal networks in the model proposed here shows that the nutrient dynamics described by the SODE trajectories (9), obtained by the TA method, now provide trajectories that model these dynamics.

As in the first SODE (1) where classical ecological interactions are modeled and the system provides trajectories that model the different dynamics under such interactions, SODE (9) described here provides trajectories that will describe the effect of carbon cycling, providing different interaction schemes according to the coordinates and assumed values for their coefficients. Although for different choices of these elements the dynamics effects can change, the geometric concepts of the TA [27] derived from this system allow the study of the stability, among other aspects, of the system according to this approach.

Using analytical trophodynamics [15, 16], a model was presented for two tree species (dominant and codominant) inhabiting a fixed region with a specific microclimate, defined by means of solar incidence, wind speed, air temperature, substrate, and pH. In this model, it was considered species interacting with the plant community in the local understory and with microorganisms from the rhizosphere, and of the subsoil of the soil-root interface [46, 47]. The results obtained extended other existing ones [20], in which the interaction coefficients had zero curvature. The new model constructed [23] took into account the potential for arbuscular infectivity defined by the mycorrhizal network, which contributes to the understanding of the interactions between trees and plants in the same region. Commonly, these interactions are complex, not being limited to the most usual in ecology: predation, parasitism, and symbiosis. Mycorrhizal networks, present in practically all subsoil, connect at least two local plants through the mycelia, which serves as a physical medium for the transmission of nutrients [1, 8]. This feature was also contemplated by the proposed mathematical model. It is also noteworthy that the results obtained corroborate experimental data [18, 44] that describe two species of trees – namely, the Douglas-fir (Pseudotsuga menziesii) and the paper birch (Betula papyrifera) – present in the same forest region of Canada’s British Columbia. Thus, it was possible to mathematically investigate the nutrient dynamics in addition to classifying possible stability scenarios [15]. Additionally, through the concept of Eco-strain [20], the forest dynamics were studied [23] under the addition of external environmental effects that influence the flux of nutrients. Conceptually, this aspect is formulated by analogy with the mechanics of the continuous medium. Including the infectivity potential in the proposed model also allowed us to understand the plastic and elastic deformations in the considered allometric production space.

We emphasize that the model proposed in [23] also allow investigating interactions in scenarios that include rhizosphere microorganisms, abiotic soil variables, grafts, and mycorrhizal network. Through algebraic computation [22, 30], experiments were carried out that contributed to the validation of the results. In the Appendix, there are codes of computational experiments related to the following application examples.

Application 1. Assume that the two species of trees in a forest are both close as Douglas-fir and paper-birch in forests of British Columbia [44]. As mentioned above, the dynamic of ecological interactions among these involve the understories and mycorrhizal components. Using the indices $\#1$ to denote Douglas-fir and $\#2$ for paper-birch consider the following system [23], derived from the Huxley space defined by $F_{\alpha }$ and $\alpha =1$, namely $F_1=\sqrt{{\left(y^1\right)}^2+{\left(y^2\right)}^2}$, that correspond to the Randers geodesics,

$$\begin{array}{l}\frac{d^2{X}^1}{{\mathit{ds}}^2}={\mathrm{\Omega }}_{12}{Y}^2,\\ \frac{d^2{X}^2}{{\mathit{ds}}^2}=-{\mathrm{\Omega }}_{12}{Y}^1,\end{array}$$E13

for $Y^i=\frac{{\mathit{dX}}^i}{\mathit{ds}}$, $b_1=-c_2\left(x^2-x_{\ast }^2\right)$, $b_2=c_1\left(x^1-x_{\ast }^1\right)$, ${\mathrm{\Omega }}_{12}=c_1+c_2$, and $c_1,c_2$ are constants. Therefore,

$${\left(Y^1\right)}^2+{\left(Y^2\right)}^2=R^2,$$E14

where $R$ is a constant. This solution shows that the net carbon flux varies periodically corroborating the biological experiments with these trees.

Application 2. Now, consider the Finsler functional defined by $\alpha =2$, $F_2=\sqrt{{\left(y^1\right)}^2+\frac{2}{\lambda }+{\left(y^2\right)}^{\frac{2}{\lambda }}}$, so that, using normal coordinates ${\overline{x}}^i$, with rates $\frac{d{\overline{x}}^i}{\mathit{dp}}=N^i$ and the uniform condition to insure that $\frac{N^i}{N^j}$ are constants over a time interval or close to climax equilibrium, the Randers geodesics are given by

$$\begin{array}{l}\frac{{\mathit{dN}}^1}{\mathit{dp}}={\mathrm{\Omega }}_{12}\left(-g^{12}{N}^1+g^{11}{N}^2\right),\\ \frac{{\mathit{dN}}^2}{\mathit{dp}}={\mathrm{\Omega }}_{12}\left(-g^{22}{N}^1+g^{12}{N}^2\right).\end{array}$$E15

In [23], it is shown that this system is also periodic both in $p$ parameter and in real-time $t$. Moreover, the $N$-constant ratio occurs in original variables, $\frac{{\mathit{dx}}^1}{{\mathit{dx}}^2}$, but with different constants.

These two applications of the model presented here illustrate how the AT method allows a mathematical description of complex dynamics such as those involving plants under different effects of the external environment. The complex mycorrhiza’s role in the interactions of plants and its effects on neighboring species, investigated under an analytical approach, can be tested from a statistical perspective in the next step.

Acknowledgments

The second author thanks Prof. Peter Antonelli and Prof. Solange Rutz for introducing him to such a rich and exciting area of research.

Conflict of interest

The authors declare no conflict of interest.

Computational experiments with FINSLER package [22] Figures 1 and 2.

Figure 1.

Application 1 experiment codes.

Figure 2.

Application 2 experiment codes.