Numerical Modeling of Chemical Compounds’ Fate and Kinetics in Living Organisms: An Inverse Numerical Method for Rate Estimation from Concentration

Abstract

Emerging chemical compounds are ubiquitous in all environmental compartments and may pose a risk to biota ecosystems. The quantification and prediction of environmental partitioning of these chemicals in various environmental compartment systems (water, sediments, soil, air, biota) is an important step in the comprehensive assessment of their sources, fates, and not finally of their uptake potential by various living organisms of ecosystems.

1. Introduction

Thousands of chemicals are used in industry, agriculture, pharmacy, commerce, and daily life. With that, a large number of chemical compounds enter the environment. Often these are considered with potential harmful effects on environmental media quality and biota safety. Thus, the understanding of them is crucial both for a better management of their use and for the better protection of the environment and living organisms. Monitoring these chemical compounds frequently is time-consuming and requires large financial efforts; one of the most cost-effective as well as time-efficient methods of evaluating their behavior in the environment and living organisms could be the use of predictive numerical models.

Experimental data and models for chemical compounds’ fate and kinetics in living organisms play a crucial role for assessing the potential human and ecological risks associated with chemical use.

Plants are receptor organisms and could be either direct or indirect vectors for chemical exposure to all other organisms [1]. In the first instance, the generated experimental data considering chemical concentrations in different media of the environment and biota are necessary to improve our understanding on plant-chemical-environment interactions. These, in turn, admit and bring forward the development of better scientific knowledge as well as conceptual and predictive models on chemical partition, fate, and uptake [2]. The strong interconnections between experimental data and model development are continuous and a long-term updated process which is needed to advance our ability to provide reliable quality information that can be used in various environmental protection contexts and regulatory risk assessments [1].

At this moment there are no standard protocols both for chemical compounds’ bioaccumulation data generation and for data use for prediction through numerical methods [1, 3]. For the reliable modeling of plant-chemical-environment interactions with the major goal to predict chemical compounds’ fates and kinetic in living organisms, it is necessary to understand and keep into account all process, phenomena, and characteristics of both chemicals and receptors (living organisms and the environment) and the interconnected process between them. Inconsistent data collection, inaccurate generation of them, or reporting them with gaps will provide improper and less useful information for their application in assessment and numerical model development.

This chapter is about to find optimal numerical modeling ways considering chemical compounds’ fate and kinetics in a living organism, specifically plants. The aim of this study was to propose a numerical procedure which estimates the highest accumulation rate of a chemical compound of interest for a growing living organism and to validate the procedure.

2. Problem formulation

1. The presence of chemical compounds in the environment and biota is the subject of global interest with the general aim to ensure that their impact on humans, other living organisms, and the environment is minimized. Chemicals that once enter the environment can persist or support different types of transformations resulting in new transformation products. Based on the characteristics of transformation products, these transformation processes could be either of benefit (attenuation processes) or be hazardous (when the resultant transformation products possess more hazardous characteristics than the “parent” compound).

2. Chemical compounds’ transformation in the environment: Based on chemical compounds’ specific physicochemical properties as well of biogeochemical and physicochemical characteristics of the media in which they are discharged (soil, water, air), these chemicals can be distributed across different environmental compartments (soil, surface and/or underground water , air, sediments, etc.) and biota [4]. In almost all cases their accumulation in environmental or biota compartments is characterized by continuous dynamic processes such as volatilization, degradation, precipitation, sorption, and so on, processes that often have the potential to end up in the formation of new chemical compounds called as transformation products of the initial chemical compound (Figure 1). Often these compounds could enter in new reaction processes and for other new end products. A schematic diagram of principal processes that could take place in both environment and biota and their interconnection with the “initial” chemical compounds and “resultant” transformation products is presented in Figure 1.

Figure 1.

A schematic presentation of the main processes involved in chemical compounds and corresponding potential transformation of products’ fate between different environmental compartments and biota.

Challenges considering potential transformation products are given by their physicochemical properties that in several cases could make them more harmful against environment quality or biota safety than the initial chemical compound. Both processes as well as the resultant transformation products’ formation are directly dependent on the environmental conditions as well properties of chemicals.

1. Chemical compounds’ uptake by plants: Chemicals’ uptake by plant organisms is a system of complex and multi-step processes. These processes could be classified firstly as chemical uptake and transportation between different anatomical compartments (e.g., root to any other anatomical compartment) and secondly as chemical uptake from different environmental compartments (route of exposure) and plant anatomical compartments (particle deposition, vapor uptake from the atmosphere, and so on). The amplitude of these processes is determined by physicochemical properties of the chemical that is under uptake [5]. Current literature presents clear-cut evidence that the availability of most organic chemical compounds is governed on the one hand by their lipophilicity and on the other hand depends on the organic matter (OM) content of the soil under consideration [6]. Some compounds form “bound” residues with organic matter (OM) or humus particles in the soil. Besides, the nature and rooting pattern of the vegetation will have greater influence on the solubility of chemicals. Exuding up to 25% of the net carbon fixed during photosynthesis into the rhizosphere, plants modify given soil-chemical interactions in multiple ways. Secondary plant products (phenolic) and soil bioactive compounds (carbohydrates, organic acids, etc.) could also impact soil micro-biodiversity that could influence in a positive way transformation of organic pollutants to reactive metabolites [7]. For example, it has been demonstrated that isoproturon is metabolized to available plant and reactive compounds in rhizosphere soil [8], while the bacterial conversion of arochlors to reactive metabolites has been one of the early results of bioremediation studies [9].

2. Probably, one of the most effective ways to study chemical behavior and fate is to use mathematical fate models. Mechanistic environmental models use mathematical equations which describe the parameters of an environment (e.g. data on flows, depths, pH, temperature, etc.) interconnected with the physicochemical properties of the chemical compounds under various conditions with the final aim of predicting their fate in the environment. According to [10], this can be an inexpensive and suitable approach for setting the limits for discharges in the environment of certain chemical compounds, and since the initial parameter description has been set up and validated by in situ and laboratory data, it can be studied with a minimum set of analysis (e.g., only the quantification of chemical compound inputs to the ecosystem) [10].

3. Framework for chemical compound interaction with the environment and living organisms

The ability of numerical models to accurately predict concentrations of target chemical compounds in any living organism depends on the model’s ability to mimic the processes involved in their uptake, and this must be assessed before they can be confidently applied [10]. After that it is necessary to consider all of these processes in order to include them in the numerical model that wants to be developed [11].

Soil-root transport: The uptake of chemicals by the root from the soil is mediated in high percentage by soil water content through the plant transpiration process [12]. A large number of organic chemicals also can be sorbed or bound to the components in soil (clay, iron oxides, organic matter), those often found in the rhizosphere in significant amounts [13, 14]. Also, lipophilic organic chemicals possess a greater tendency to partition into plant roots than hydrophilic chemicals. Although chemicals properties are important predictors of the uptake potential, the physiology and composition of the plant root itself is also a significant influence, with differences in the uptake potential explained by the varying types and amounts of lipids in root cells [15]. Uptake from the external media is often expressed as a root concentration factor (RCF), which is the ratio of chemical concentration in the root to the concentration found in external media [13, 16].

Transfer from roots to other anatomical compartments of plants: The major factor that illustrates the amount of a chemical compound that was transferred from the plant root part to other anatomical compartments is the transpiration stream concentration factor (TSCF) which is the ratio of chemical concentration in transpiration to the concentration found in the external part. TSCF could be predicted from knowledge of the chemical compound lipophilicity, with maximum uptake, a logK_OW about 1 [17]. Once the chemical is transported to the stem, plant water and solutes take it and continue to transport these chemicals to the rest of the anatomical compartments through vascular systems and cell tissues [18].

Vapor or gas uptake from ambient air: Another exposure route with chemicals in case of plant materials could be the ambient air that contains a large number of contaminants. This exposure route is governed by gaseous exchange and facilitates the transport and uptake of chemicals that are volatiles and which are more easily partitioned in air than in water. This has been shown to be the main uptake pathway in the above-ground plant parts for a variety of chemical compounds (e.g. PCBs, tetra- and hexa-chlorinated PCDD/Fs) [19]. Previous studies have reported a good correlation between shoot uptake and chemical properties of compounds (K_OW, Henry’s Law constant, octanol-air partition coefficient, etc.). In studies presented by [20], it was evidenced that gaseous uptake is the primary pathway for chemicals with an octanol-air partition coefficient (logK_OA) less than 11 [20].

Particulate deposition on plant surfaces: Pesticides as well as other chemical contaminants are bound to soil particles which may be transported by wind and/or rain and deposited on the above-ground anatomical compartments of plants. Studies presented by [5, 21], evidenced that dry deposition onto the leaf of suspended particles that contain PCCD/Fs is the major route of uptake due to PCCD/Fs permeation through the cuticle. Similarly, in their studies, wet deposition was shown which could also be the dominant deposition mechanism for organic chemicals with Henry’s Law constant of less than $1\times {10}^{-6}$ [5, 21].

To date, a number of mathematical models have been developed to facilitate the exposure assessment of chemical contaminants, with important results in the modeling of pollutants’ multimedia fate and the modeling of pollutants’ linkage with transformation products, especially in water environmental compartments [4].

4. Case study presentation

The properties of wild growing mushrooms make them valuable resources both in culinary practices and in pharmaceutical practices. They are recognized as healthy food with low contents of calories and fats but high in vitamins, minerals, and vegetable proteins. Their suitability for use by the pharmaceutical industry is given by their rich antioxidant chemical constituents that are capable of preventing the human body from oxidative damage [22]. It is also known that mushrooms could be considered as good bioindicators for the evaluation of environmental pollution, since they are known to accumulate a broad range of chemical compounds [23]. The aim of this study was to propose a numerical procedure which estimates the highest accumulation rate (R) of a chemical compound on the entire anatomical compartments of a mushroom body. Such data could lead to improvement in both food quality assurance and environment safety assessment. Analytical assessments on mushroom samples have shown that the accumulation potential of chemical compounds varies with mushroom species and varieties and also varies between the same mushroom anatomical compartments (see Figure 2) as well between mushroom development stages (“age”).

Figure 2.

Chemical compound concentration variations in different anatomical compartments: (a) concentration variation in the first anatomical compartment (basal bulb) of the mushroom; (b) concentration variation in the second anatomical compartment (stipe) of the mushroom; (c) concentration variation in the third anatomical compartment (cap) of the mushroom.

5. Solving modality path selection and motivation: the inverse numerical method for rate estimation from concentration

Inverse problems are extremely frequent in interdisciplinary science subjects. A large scale of mathematical and numerical techniques for solving scattering problems as well as other inverse problems usually exist [24]. These methods are often very different from the methods used for solving direct problems due to the differences in mathematical structure and input data [25].

In our study, the estimation process of the chemical compound accumulation rate was built on the following differential equation:

$$\frac{d\left(\mathit{p\varphi }\left(z\right)\mathit{\chi C}\left(z\right)\right)}{\mathit{dz}}-\frac{d}{\mathit{dz}}\left[\mathit{p\varphi }\left(z\right)\left(W_C-\overline{\mathit{ET}}\right)\frac{\mathit{dC}\left(z\right)}{\mathit{dz}}\right]+\mathit{p\varphi }\left(z\right){\mathit{Up}}_f\left(C\left(z\right)-{\mathit{Ba}}_f\right)=R$$E1

where the elements which may affect the rates are given as follows: $z$ is the height [mm], $\mathrm{\Phi }$ is the diameter [mm], $p$ is the porosity, $\chi$ is the hydration coefficient, $W_c$ the is saturated hydration factor, $\overline{\mathit{ET}}$ is the evaporation-transpiration coefficient, ${\mathit{Up}}_f$ is the uptake factor, ${\mathit{Ba}}_f$ is the bioaccumulation factor, $R$ is the rate of accumulation for the chemical compound of interest, $C$ is the concentration data of the target chemical compound [ng⋅g⁻¹].

The rate estimation model had as a starting point the one-dimensional transport-reaction equation for dissolved compounds presented by Lettmann et al. [26]. In this chapter the same type of equation was used but this time the equation is based on the main factors that can influence in some way the assimilation rate of a chemical compound in a vegetal organism—specifically in a mushroom body. Our model was supported by concentration data (C) of the target compounds, which were measured in the laboratory from cross-sections taken at every 2 mm over the whole body of the studied mushroom species. Also, the concentration measurements correspond to cross-sections taken at every 2 mm, and each section was divided into three concentric subintervals with regard to diameter.

The goal of the first step is to approximate $R$ using the left-hand side of Eq. (1). The approximation of differential operators from the left side of the proposed equation has been solved using smoothing spline functions [27, 28].

Model validation was performed by solving a two-point boundary differential equation relative to Eq. (1) on the interval given by extreme values of z and comparing with the measured values of C. The numerical method for target compounds’ accumulation rate validation (the solution of BVP) was implemented using MATLAB¹ bvp4c function [29, 30]. Good concordance was identified between the measured concentration and the concentrations computed by the solution of Eq. (1), given the rate R estimated during the validation process. The concordance is given by mean-square deviation. In the paper presented by Lettmann et al., [26] the approximation of the linear differential operator is performed by finite differences, while in our case its approximation was done through smoothing spline. Also, the rate estimation was generated randomly while in our case rate estimation was based on experimental data obtained in the laboratory. Their work, due to the nature of the practical problem, has no constraint on volume while in our case we were limited to the relative small volume and dimensions of the studied mushroom species.

6. Analysis and modeling: the estimation of R

6.1. Input data

We select for the study mushrooms from species Macrolepiota procera, which is one of the most popular in consumption and frequency from our country. Experimental measurements of major parameters involved in our model (parameters from Eq. (1)) were done on mature samples collected from the natural habitat where the evidence of potential contamination with chemical compounds exist. The measured data were grouped by compartments—basal bulb, stipe, gills and cap—representing the main anatomical parts of the mushroom body. Their form is illustrated in Table 1.

Table 1.

Input parameters and input data structure.

To convert our problem from bi-dimensional to a one-dimensional one, we consider a weighted mean for concentrations and parameters. The piecewise constant parameters ($p$, $\chi$, $W_c$, $\overline{\mathit{ET}}$, ${\mathit{Up}}_f$, ${\mathit{Ba}}_f$) are weighted by height (z), and the concentrations are weighted by diameter, on sections approximately orthogonal to median axes (see Figure 3).

Figure 3.

Schematic representation of target chemical compound concentration distribution around mushroom anatomical compartment.

6.2. Differential operator and its approximation

Once we have the averaged concentration we compute the rate using the formula:

$$R\left(z\right)=\frac{d\left(p\mathrm{\Phi }\left(z\right)\mathit{\chi C}\left(z\right)\right)}{\mathit{dz}}-\frac{d}{\mathit{dz}}\left[p\mathrm{\Phi }\left(z\right)\left(W_c-\overline{\mathit{ET}}\right)\frac{\mathit{dC}\left(z\right)}{\mathit{dz}}\right]+p\mathrm{\Phi }\left(\mathrm{z}\right){\mathrm{Up}}_{\mathrm{f}}\left(\mathrm{C}\left(\mathrm{z}\right)-{\mathrm{BA}}_{\mathrm{f}}\right)$$E2

The numerical differentiations involved in Eq. (2) are critical operations, leading to large errors. Lettmann et al. [26] performed using finite differences methods followed by a Tikhonov least-squares regularization [31, 32, 33]. Our approach is different and is based on smoothing splines. The diameter $\mathrm{\Phi }$ is approximated by a cubic piecewise Hermite spline (MATLAB function pchip) and the concentration by a smoothing spline (MATLAB function spaps, in the spline toolbox or in the curve fitting toolbox in newer versions). Since our approximations are piecewise polynomial, the computation of their derivatives is straightforward (using fnval and fnder functions) [34, 35]. The utilization of the smoothing spline for concentration allows us to reduce the propagated errors and to perform a correction equivalent to Tikhonov regularization [34].

6.3. The smoothing spline

We look for a spline function f, in the B-spline basis, that minimizes the expression: $\mathrm{\rho }{\mathrm{E}}_{\mathrm{f}}+\mathrm{F}\left({\mathrm{D}}^{\mathrm{m}}\mathrm{f}\right),$ where $E_f$ is the distance of the spline function f from the given data, given by:

$$E_f={\displaystyle \sum _{j=1}^n}{w}_j{‖y_j-f\left(x_j\right)‖}^2,$$E3

$\mathrm{F}\left({\mathrm{D}}^{\mathrm{m}}\mathrm{f}\right)$ is:

$$F\left(D^{m}f\right)=\underset{x_{\mathit{min}}}{\overset{x_{\mathit{max}}}{\int }}\lambda \left(t\right){‖D^{m}f\left(t\right)‖}^2\mathit{dt},$$E4

and ρ is the smoothing factor. The spaps function uses the algorithm described in Reinsch’s work [32]. For additional details on smoothing and interpolation splines, see deBoor’s book [33] and the MATLAB curve fitting toolbox user’s guide [36, 37].

7. Validation and solution of the differential equation

This step has a double purpose: to prove that the approximation of rate, given by Eq. (2), is sufficiently accurate and to compute the concentration from the rate, without performing any measurement. We want to solve the two-point boundary value problems Eq. (1) and boundary conditions:

$$\left\{\begin{array}{c}C\left(z_{\mathit{\text{initial}}}\right)=C_{0,}\\ C\left(z_{\mathit{\text{final}}}\right)=C_n.\end{array}\right.$$E5

Our solution uses the collocation method [38]. The independent variable z in Eq. (1) means the length of the path along the medial axis.

8. Discussions and main conclusions

Physiological events modeling, as uptake, bioaccumulation, or metabolism, and so on, in living organisms are extremely difficult both due to the complex nature of physiological processes and due to the complexity of the biological system that is modeled. For this reason the existence of implemented models in this area is very scarce, if not almost inexistent. Most existing models from the literature refer to models which are applied to a micro-scale fragment from a biological system (intercellular models) and less for globalized macro-scales that integrate multiple events [29, 39].

In this chapter we tried to overcome this challenge, trying to model the accumulation rate for a chemical compound based on experimental data obtained in the laboratory after analysis conducted on numerous mushroom specimen analyses of Macrolepiota procera. Our first approach was on the anatomical compartments of the mushroom body, but the results were not satisfactory.

Initially, the independent variable z was the height of the same compartment. Because compartments as cap and lamellae had an insufficient number of data (since the length of cap and lamellae is 1.2 cm and the minimal width for sample collection was from the section taken from 2 and 2 mm), we obtained large deviations between measured concentrations and computed concentrations at inter-compartment boundaries (see Figure 4).

Figure 4.

The averaging of parameters (ODE solution, smoothing spline) and target chemical concentrations on the whole mushroom.

Thus, because in chemistry in a given volume the concentration of a chemical compounds is the same in any point of these volume, we considered in the next that we have several points of concentration data of same value in horizontal sections of the cap (Figure 3).

These drawbacks lead us to modify the approach mentioned in Section 6 on Analysis and modeling: The estimation of R. Due to symmetry we considered a half of an axial section and a medial axis of the section. Now z is the length of the path on medial axis. To apply the classical theory on ordinary differential equations to Eq. (1), we need to have a function of class $C^2$ on the domain of z, while the parameters defining the rate are piecewise constant. For this reason we considered the weighted average of these parameters on the whole length of the medial axis. The diameter $\Phi \left(z\right)$ was approximated by the piecewise Hermite cubic spline of measured diameters. Our choice is motivated by the fact that these splines are shape-preserving. The next step was the application of Eq. (2) to compute the rate R. The rate R, computed along the medial axis, as described in the previous section, is plotted in Figure 5.

Figure 5.

The graph of rate R.

Analyzing data obtained for monitored chemical compound rate, it was possible to observe that larger fluctuations are present in the mushroom stipe while in the caps part (cap and lamellae) a decreasing tendency is registered. These data are in correlation both with infield experimental measurements and with the computed concentration obtained from our model—see Figure 6 where the concentration C is presented, after the solution of the differential equation (1) with boundary conditions (5).

Figure 6.

The graph of concentration C, obtained by the solution of a two-point boundary value problem.

In order to assess the accuracy of our model we plot the initial and the computed concentration data (Figure 7).

Figure 7.

The graph of initial and computed concentrations.

There is good matching, as Figure 8 shows. The least square deviation for the concentration is 5.1384e-005, and it was computed at points corresponding to the measured z and concentrations.

Figure 8.

Differences between the initial and the computed concentration.

The model proposed for the rate estimation of a chemical compounds in living organisms, specifically a mushroom, is new. Based on our knowledge, up to date, there is no paper on the rate estimation for a specific chemical compound in a vegetal system. The model presented by Lettmann et al. [26] has completely different premises—there is no volume limitation in their case study while in our case we were limited to the smaller dimension of the studied living organism, the mushroom species Macrolepiota procera. Once have the rate for a species, we can compute concentrations via the solution of Eqs. (1) + (5) without doing any additional laboratory measurement.

Acknowledgments

This chapter was funded by the Core Program, under the support of ANCS, project OPTRONICA V number PN PN18-28 02 0.