Fifth Mississippi State Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, Conference 10, 2003, pp 89-99. http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu (login: ftp)

Modelling the interactions of mixtures of organophosphorus insecticides with cholinesterase^{∗}

Russell. L. Carr, Howard W. Chambers, Janice E. Chambers, Seth F. Oppenheimer, & Jason R. Richardson

Abstract

The organophosphorus (OP) insecticides are one of the most widely used and important insecticide classes. These insecticides exert toxicity through inhibition of the critical nervous system enzyme cholinesterase (ChE) which functions to rapidly destroy the ubiquitous neurotransmitter acetylcholine. When ChE is inhibited, the acetylcholine accumulates, causing hyperactivity within the cholinergic pathways. Considerable effort has gone into assessing the risks of various OP insecticides. Unfortunately, people are often exposed to different OP insecticides in different dosages at different or overlapping times. The usual statistical methods seem inadequate to the task of assessing the effect of OP mixtures. This paper will discuss a simple model using systems of ordinary differential equations. Using this model, we have had success in predicting the effect of cumulative in vitro OP compound exposure in terms of ChE inhibition using data from experiments measuring ChE inhibition by a single OP compound. We will describe our model and compare our simulations to in vitro experiments where binary mixtures have been used.

1 Introduction

The organophosphorus (OP) insecticides are one of the most widely used and important insecticide classes. These insecticides exert toxicity through inhibition of the critical nervous system enzyme cholinesterase (ChE) which functions to rapidly destroy the ubiquitous neurotransmitter acetylcholine. Inhibition of ChE by OP compounds is through covalent bond formation, and the inhibited ChE is quite persistent (half life of reactivation to uninhibited ChE is hours to days). An additional covalent reaction can occur, termed aging, which renders the ChE molecule permanently inhibited and incapable of recovering enzymatic activity. When ChE is inhibited, the acetylcholine accumulates, causing hyperactivity within the cholinergic pathways. Considerable effort has gone into assessing the risks of various OP insecticides. Unfortunately, people are often exposed to different OP insecticides in different dosages at different or overlapping times. The usual statistical methods seem inadequate to the task of assessing the effect of OP mixtures [6, 7, 8, 9, 10, 11, 12, 16]. This paper will discuss a simple model using systems of ordinary differential equations. Our approach will be somewhat different than those used previously in the toxicology literature [3]. Using this model, we have had success in predicting the effect of cumulative in vitro OP compound exposure in terms of ChE inhibition using data from experiments measuring ChE inhibition by a single OP compound.

∗Mathematics Subject Classifications: 92C45.

Key words: Reaction kinetics, reaction modelling \begin{matrix}c \\ ◯ \end{matrix}2003 Southwest Texas State University.

Published February 28, 2003.

This work is supported by a grant from the American Chemistry Council.

In the second section of this work, we will describe our model. Following the model description, the next section will give the calibration results. In the fourth section, we will compare our simulations to in vitro experiments of simultaneous exposures to two OP inhibitors. The agreement seen in Section Four is excellent. The fifth section briefly discusses time asymptotic results. Full information on the experimental methods and results are forthcoming [1].

2 The model

We will start by discussing the usual modelling approach.

The usual models for chemical kinetics involve elementary reactions of the following form [2, 4]:

aA + bB → cC + dD,

where the uppercase letters represent concentrations in moles per liter or a similar set of units and the lower case letters are natural numbers. In this case of an elementary reaction, the standard model is given by the following differential equations:

\frac{dA}{dt} = −akA^{a}B^{b}
\frac{dB}{dt} = −bkA^{a}B^{b}
\frac{dC}{dt} = ckA^{a}B^{b}
\frac{dD}{dt} = dkA^{a}B^{b}.

That is, the rate of reaction is always proportional to products of the concentration raised to the number of molecules that participate in each reaction. The overall order of the reaction, a+b, is also the molecularity of the reaction, where molecularity is number of molecules taking part in the reaction. It has been observed by some investigators that not all such reactions are quite so simple. It may be because of intermediate reactions or multiple reaction paths, but the above model can fail. An alternative model is proposed in [14] of the form

\frac{dA}{dt} = −akA^{α}B^{β}

EJDE/Conf/10 Carr, Chambers, Chambers, Oppenheimer, & Richardson 91

\frac{dB}{dt} = −bkA^{α}B^{β}
\frac{dC}{dt} = ckA^{α}B^{β}
\frac{dD}{dt} = dkA^{α}B^{β}

where the exponents, called the partial orders, do not necessarily have a relationship with the molecularity and are in fact obtained empirically. This is the modeling approach we will follow. We will consider single reactions first. Let us consider our case of ChE and a variety of inhibitors. Let c be the molarity in solution of ChE at time t and x_{i} be the molarity of the i^{th} inhibitor at time t. We will model the reaction of ChE with the inhibitor by the following system of differential equations

\frac{dc}{dt} = −k_{i}c^{α_{i}}x_{i}^{β_{i}}
\frac{dx_{i}}{dt} = −k_{i}c^{α_{i}}x_{i}^{β_{i}}
c(0) = c^{0}
x_{i}(0) = x_{i}^{0}

where k_{i} is called the rate coefficient and α_{i} and β_{i} are the partial orders. These constants are all empirically derived from the single inhibitor experiments. The reader will observe that we are not modelling the concentration of the result of the reaction. This is because the rate of the back reaction (reactivation) is so slow compared to the time scale of the experiments, fifteen to thirty minutes, that the back reaction will have negligible effect on the experimental results. We will denote the solution of equation (1) as c (α_{i}, β_{i}, k_{i}, x_{i}^{0}; t) .

The experimental data used in this work are in the form of inhibition curves. That is, the experimentalists will start with several samples containing a fixed concentration of ChE. Different amounts of an inhibitor are added and the fraction of the ChE which is inhibited is measured at fifteen minutes. That is, the data sets consist of several different initial concentrations of inhibitor i, {x_{i1}^{0}, …, x_{iN}^{0}} and the percentage of ChE inhibited after fifteen minutes I_{ij}. Here, x_{ij}^{0} is the initial concentration of inhibitor i in the j^{th} experiment. For each inhibitor, fifteen experiments were carried out and N = 15. (We observe that this is quite a simplification of the actual experimental procedure and the reader is referred to [1, 13].) The method of determining the unknown constants α_{i}, β_{i}, and k_{i} was as follows. We used the built–in numerical program ODE45 in the Mathlab programming language [15] to approximate the solution of equation (1) with initial inhibitor concentration x_{ij}^{0} at time 15 minutes, c(α_{i}, β_{i}, k_{i}, x_{ij}^{0}; 15), and compute the percent inhibition at 15 minutes for this set of parameters

p (α_{i}, β_{i}, k_{i}, x_{ij}^{0}; 15) = (1 − \frac{c(α_{i}, β_{i}, k_{i}, x_{ij}^{0}; 15)}{c^{0}})100.

For inhibitor i we define the objective function

O (α_{i}, β_{i}, k_{i}) = ∑_{j=1}^{15} (I_{ij} − p (α_{i}, β_{i}, k_{i}, x_{ij}^{0}; 15))^{2} .

We find the parameters α_{i}, β_{i}, and k_{i} by minimizing the objective function using the Matlab function fmins.

Once the constants have been found, we can give the model for any combination of two inhibitors:

\frac{dc}{dt} = −k_{i}c^{α_{i}}x_{i}^{β_{i}} − k_{n}c^{α_{n}}x_{n}^{β_{n}}  (1)
\frac{dx_{i}}{dt} = −k_{i}c^{α_{i}}x_{i}^{β_{i}}
\frac{dx_{n}}{dt} = −k_{n}c^{α_{n}}x_{n}^{β_{n}}
c(0) = c^{0}
x_{i}(0) = x_{i}^{0}
x_{n}(0) = x_{n}^{0}

This is the model we will use below for the simultaneous exposures. We note that this model is much more general than we indicate here. Using it in slightly different form, we can handle any number of inhibitors, any exposure regime, reactivation, and aging. If the various inhibitors react with each other or other more complex interactions occur, these basic models still form the building blocks of the model.

We observe that in the work below the initial ChE concentration was not measured directly, but was obtained using data in [13] and a published abstract [5]. The value used is c_{0} = 4.045 × 10^{−13} moles per liter.

3 Model Calibration

We will report on the calibration of the model for two OP inhibitors in this section, chlorpyrifos-oxon and paraoxon. As noted in the previous section, the unknown parameters are found by minimizing the objective function O (α_{i}, β_{i}, k_{i}). We will subscript the parameters associated with chlorpyrifos-oxon with C and we will subscript the parameters associated with paraoxon with P.

Below, we will give the identified parameters, a chart comparing the inhibition predicted by the calibrated model, and a graph of the model inhibition curve with the experimental data points.

Chlorpyrifos-oxon

The identified parameters for chlorpyrifos-oxon are

k_{C} = .0498

α_{C} = 1.1616
β_{C} = .9670

We now give a chart comparing the model’s predictions with the experimental data. Observe that the experimental error can be on the order of ±3% inhibition:

Nanomoles of chlorpyrifosoxon per liter at time 0 Percent inhibition of ChE observed at 15 minutes Percent inhibition of ChE predicted at 15 minutes Predicted − observed

.5  8.7  10.2  -1.5
.5  9.0  10.2  -1.2
.5  9.6  10.2  -.6
.5  9.97  10.2  -.23
.5  10.3  10.2  .1
1  17.7  18.8  -1.1
1  19.6  18.8  .8
1  19.7  18.8  .9
1  20.2  18.8  1.4
1  20.9  18.8  2.1
1.7  27.0  29.1  -2.1
1.7  28.0  29.1  -1.1
1.7  28.5  29.1  -.6
1.7  28.8  29.1  -.3
1.7  31.3  29.1  2.2

The information may be summarized statistically as follows. The maximum of the absolute error is 2.2, the mean error is −.082, the median error is −.3, and the standard deviation in the error is 1.1308. All of the units are percent inhibition.

The graph of the predicted fifteen minute inhibition curve along with the experimental data may be seen in Figure 1.

Paraoxon

The identified parameters for paraoxon are

k_{P} = .0050
α_{P} = 1.1160
β_{P} = 1.0133

We now give a chart comparing the model’s predictions with the experimental data. Observe that the experimental error can be on the order of ±3% inhibition:

Nanomoles of paraoxon per liter at time 0 Percent inhibition of ChE observed at 15 minutes Percent inhibition of ChE predicted at 15 minutes Predicted − observed

3.5  9.7  10.16  0.46
3.5  9.7  10.16  0.46
3.5  9.97  10.16  0.19
3.5  10.0  10.16  0.16
3.5  11.7  10.16  -1.54
7.5  19.6  20.57  0.97
7.5  20.0  20.57  0.57
7.5  20.8  20.57  -0.23
7.5  21.0  20.57  -0.43
7.5  22.0  20.57  -1.43
12  28.8  30.78  1.98
12  30.5  30.78  0.28
12  31.0  30.78  -0.22
12  31.3  30.78  -0.52
12  33  30.78  -2.22

The information may be summarized statistically as follows. The maximum of the absolute error is 2.22, the mean error is 0.1013, the median error is −0.16, and the standard deviation in the error is 1.0521. All of the units are percent inhibition.

The graph of the predicted fifteen minute inhibition curve along with the experimental data may be seen in Figure 2.

4 Simultaneous exposure to a binary mixture

Using the parameters obtained above, we will use the model

\frac{dc}{dt} = −k_{P} c^{α_{P}} x_{P}^{β_{P}} − k_{C}c^{α_{C}} x_{n}^{β_{C}}  (2)
\frac{dx_{P}}{dt} = −k_{P} c^{α_{P}} x_{P}^{β_{P}}
\frac{dx_{C}}{dt} = −k_{C}c^{α_{C}} x_{n}^{β_{C}}
c(0) = c^{0}
x_{P} (0) = x_{P}^{0}
x_{C}(0) = x_{C}^{0}

to predict the percent ChE inhibition when two inhibitors, chlorpyrifos-oxon and paraoxon, are present.

Using this model we got excellent agreement with experiment. By excellent, we mean that the inhibition predicted by the model was within experimental error of the observed inhibition. This is particularly significant as the model was calibrated independently of any of the binary mixture data. We present a table summarizing our results.

Nanomoles per liter of paraoxon at time 0 Nanomoles per liter of chlorpyrifosoxon at time 0 Observed inhibition Predicted inhibition

Observed - predicted

3.5  .5  18.1  19.17  -1.07
3.5  1  27.3  26.8  .5
3.5  1.7  36.33  35.95  .38
7.5  .5  27.1  28.41  -1.31
7.5  1  35  35.06  -.06
7.5  1.7  42  43.06  -1.06
12  .5  35  37.48  -2.48
12  1  41.7  43.19  -1.49
12  1.7  49  50.07  -1.07

The information may be summarized statistically as follows. The maximum of the absolute error is 2.48, the mean error is −0.8511, the median error is −1.07, and the standard deviation in the error is 0.9604. All of the units are percent inhibition.

The graph of the predicted fifteen minute inhibition surface along with the experimental data may be seen in Figure 3.

5 Results for larger times

Our goal in this work has been to develop models for Cholinesterase inhibition by mixtures of OP insecticides that are more accurate and flexible than those currently available. We did so in the context of standard experimental protocols for obtaining fifteen minute inhibition curves. However, we shall consider some results, both theoretical and experimental for longer times. We note that as time increases, we expect other mechanisms, such as reactivation, to come into play, thus examining this model in isolation over longer times has limited utility.

We observe for the model with only one OP

\frac{dc}{dt} = −kc^{α}x^{β}
\frac{dx_{i}}{dt} = −kc^{α}x^{β}
c(0) = c^{0}
x(0) = x^{0}

that for each t > 0, x(t) − x^{0} = c(t) − c^{0}. We may therefore write

\frac{dc}{dt} = −kc^{α} (c(t) + x^{0} − c^{0})^{β} .

If we assume that x^{0} > c^{0} as we had in all of our experiments we obtain

\frac{dc}{dt} ≤ −kc^{α} (x^{0} − c^{0})^{β} .

From this it is easy to see that if a < 1 then c reaches zero in finite time and if α ≥ 0 then lim_{t→∞} c (t) = 0. Of course, a residual of the OP will be left of the amount x^{0} − c^{0}. For the case of multiple OP’s we can construct a similar upper bound and obtain analogous results. We note that the OP for which α is smaller will dominate the reaction. If the α’s are equal, then the term with the largest k will dominate.

For longer time periods, experimental measurements become more difficult. However, we do have a comparison of the thirty minute inhibition curve for paraoxon as predicted by the model and experimental data. We will plot the model inhibition curve along with forty five data points.

The information may be summarized statistically as follows. The maximum of the absolute error is 7.6, the mean error is 2.27, the median error is 1.98, and the standard deviation in the error is 2.69. All of the units are percent inhibition. Considering the range of experimental values, this is not too bad.

Conclusion We developed a new model for enzyme inhibition by an OP insecticide. We showed that this model accurately modelled single inhibitor experiments. More significantly, using parameters obtained independently of any data from experiments using binary mixtures, the model successfully predicted the results of in vitro experiments using two inhibitors.

References

[1] J. E. Chambers, R. L. Carr, J. R. Richardson, H. W. Chambers, and S. F. Oppenheimer. The In Vitro Inhibition of Rat Brain Cholinesterase Following Simultaneous or Sequential Exposure to Binary Mixtures of Organophosphate Compounds. In preparation.

[2] I. Glassman. Combustion, Second Edition. Academic Press, Orlando, 1987.

[3] Stacey A. Kardos and Lester G. Sultatos. Interactions of the Organophosphates Paraoxon and Methyl Paraoxon with Mouse Brain Acetylcholinesterase, Toxicological Sciences Vol. 55, 118–126 (2000).

[4] James Keener and James Sneyd. Mathematical Physiology. Springer, New York, 1998.

[5] A. K. Kousba and L. G. Sultatos. Estimation of Acetylcholinesterase Concentrations for Pharmacodynamic Models. The Toxicologist Vol. 60, 241 (2001).

[6] J. P. Michaud, A. J. Gandolfi, and K. Brendel. Toxic Responses to Defined Chemical Mixtures: Mathematical Models and Experimental Designs. Life Sciences, Vol. 55 (9), 635–651 (1994).

[7] John L. Plumber and Tim G. Short. Statistical Modeling of the Effects of Drug Combinations. Journal of Pharmacological Methods Vol. 23, 297–309 (1990).

[8] Gerald Pöch, Peter Dittrich, and Sigrid Holzmann. Evaluation of Combined Effects in Dose-Response Studies by Statistical Comparison with Additive and Independent Reactions. Journal of Pharmacological Methods Vol. 24, 311–325 (1990).

[9] Gerald Pöch, Peter Dittrich, R. J. Reiffenstein, W. Lenk, and A. Schuster. Evaluation of Experimental Combined Toxicity by Use of Dose-Frequency Curves: Comparison with Theoretical Additivity as Well as Independence. Can. J. Physiol. Pharmacol. Vol. 68, 1338–1345 (1990).

[10] Gerald Pöch and Sonja N. Pancheva. Calculating Slope and ED_{50} of Additive Dose-Response Curves, and Application for These Tabulated Parameter Values. Journal of Pharmacological and Toxicological Methods Vol. 33, 137–145 (1995).

[11] Gerald Pöch, R. J. Reiffenstein, and H.-P. Baer. Quantitative Estimation of Potentiation and Antagonism by Dose Ratios Corrected for Slopes of Dose-Response Curves Deviating From One. Journal of Pharmacological and Toxicological Methods Vol. 33, 197–204 (1995).

[12] Ross L. Prentice. A Generalization of the Probit and Logit Methods for Dose Response Curves. Biometrics Vol. 32, 761-768 (1976).

[13] Jason R. Richardson, Howard W. Chambers, and Janice E. Chambers. Analysis of the Additivity of in Vitro Inhibition of Cholinesterase by Mixtures of Chlorpyrifos-oxon and Azinphos-methyl-oxon. Toxicology and Applied Pharmacology Vol. 172, 128–139 (2001).

[14] Garrison Sposito. Chemical Equilibria and Kinetics in Soils. Oxford University Press, Oxford, 1994.

[15] The Math Works Inc., The Student Edition of MATLAB, Prentice Hall, New Jersey, 1992.

[16] Aage VØlund. Application of the Four-Parameter Logistic Model to Bioassay Comparison with Slope Ratio and Parallel Line Models. Biometrics Vol. 34, 357–365 (1978).

Russell L. Carr (e-mail: rlcarr@cvm.msstate.edu) Janice E. Chambers (e-mail: chambers@cvm.msstate.edu) Center for Environmental Health Sciences, P.O. Box 6100 Mississippi State University, MS 39762, USA Howard W. Chambers Department of Entomology and Plant Pathology P.O. Box 9775, Mississippi State University, MS 39762, USA Seth F. Oppenheimer Department of Mathematics and Statistics, Drawer MA Mississippi State University, MS 39762, USA e-mail: seth@math.msstate.edu Jason R. Richardson Center for Neurodegenerative Disease, Emory University Whitehead Biomedical Research Building 575 615 Michael Street, Atlanta, GA 30322, USA e-mail: jricha3@emory.edu