Skip to main content

Potential new therapeutic modality revealed through agent-based modeling of the neuromuscular junction and acetylcholinesterase inhibition



One of the leading causes of death and illness within the agriculture industry is through unintentionally ingesting or inhaling organophosphate pesticides. OP intoxication directly inhibits acetylcholinesterase, resulting in an excitatory signaling cascade leading to fasciculation, loss of control of bodily fluids, and seizures.


Our model was developed using a discrete, rules-based modeling approach in NetLogo. This model includes acetylcholinesterase, the nicotinic acetylcholine receptor responsible for signal transduction, a single release of acetylcholine, organophosphate inhibitors, and a theoretical novel medical countermeasure. We have parameterized the system considering the molecular reaction rate constants in an agent-based approach, as opposed to apparent macroscopic rates used in differential equation models.


Our model demonstrates how the cholinergic crisis can be mitigated by therapeutic intervention with an acetylcholinesterase activator. Our model predicts signal rise rates and half-lives consistent with in vitro and in vivo data in the absence and presence of inhibitors. It also predicts the efficacy of theoretical countermeasures acting through three mechanisms: increasing catalytic turnover of acetylcholine, increasing acetylcholine binding affinity to the enzyme, and decreasing binding rates of inhibitors.


We present a model of the neuromuscular junction confirming observed acetylcholine signaling data and suggesting that developing a countermeasure capable of reducing inhibitor binding, and not activator concentration, is the most important parameter for reducing organophosphate (OP) intoxication.


Inadvertent or intentional ingestion of organophosphate (OP) pesticide is a common occurrence in agricultural areas [1] and OP nerve agents remain a threat in chemical terrorism [2]. Efforts to develop new therapeutic treatments for OP poisonings are frequently resulting in new oxime-based reactivators [3], stoichiometric bioscavengers [4], or catalytic scavengers [5]. These treatments rely upon knowledge and identification of an acute exposure, must be administered within a narrow therapeutic window, and are not broad treatments but are selective for specific OPs. Despite these limitations, the currently deployed OP therapeutics used in the U.S. (the Mark 1 pralidoxime/atropine autoinjector and pyridostigmine bromide) have been in use since the 1950s.

Current advanced computational models of OP intoxication are constructed as physiologically-based pharmacokinetic (PBPK) models in order to estimate target tissue dosimetry. These models are then connected to pharmacodynamic (PD) models to predict the biological response at the target site. Such models exist for paraoxon [6, 7], diazinon [8, 9], chlorpyrifos [9], diisopropylfluorophosphate (DFP) [7, 10], and more recently for soman [11]. The model for paraoxon poisoning was developed for the rainbow trout while the models for diazinon and DFP, in contrast, were validated against rat and human in vivo data. More broadly applicable models were developed for soman [11] and for dermal absorption of pesticides such as parathion and fenthion [12]. The primary advantage of these PBPK models is that they can provide an accurate estimate of population behaviors and predict systemic outcomes.

The work presented here develops a PD model of the mammalian neuromuscular junction (NMJ) based on an agent-based model (ABM) describing acetylcholine signaling through nicotinic receptors (Figure 1). Agent-based modeling is a discrete, rules-based method of computational modeling that focuses on the individual components of an experimental system to perform in silico experiments [13]. ABMs are well suited for cases where the modeling goal is to test the validity of a mechanistic hypothesis [14], such as the case herein where allosteric activation of AChE is proposed to protect against OP intoxication. For example, the use of an ABM to model signaling events in the NF-κB pathway showed strong correlation between ABM, differential equation approaches (ODE), and in vitro measurements [15]. Lipniacki [16] showed that a purely ODE approach within the NF-κB system does not fully account for singular events, which required superposition of stochastic modeling onto the ODE. Furthermore, a recent modeling method has extended the ABM to include even finer resolution of physical space in chemical reactions, generating a spatial model of toll-like receptor 4 immune signaling that qualitatively reproduced the observed dynamics of tumor necrosis factor secretion [17].

Figure 1
figure 1

Conceptual rendering of the neuromuscular junction and NetLogo rendering of the junction as relates to the model.

Compared with ODEs, ABM constructs are readily adapted to spatial dimensions [18]; are stochastic by nature; can easily incorporate new information by adding more agent-types or modifying rules without rewriting the entire simulation; and reproduce emergent behaviors through parallelism and stochasticity [14]. Models in the ABM paradigm can be assembled even when complete knowledge of the system being simulated is lacking, as, for example, in the case herein where the characteristics of an enzyme activator are theoretical. Finally, ABMs describe the behavior of individuals such that the simulation does not always follow the average behavior that the ODE description would provide, thus taking into account the often significant impact of “outlier events” on the overall biological process. Although the system outcome from each ABM run is different, multiple runs provide a non-parametric means to explore the variability of outcomes, including the impact of rare events, eventually converging with the ODE-based results.

Traditional ODE models can be successfully employed to predict macroscopic effects that are changing in a continual manner; however they fall short in modeling dynamic processes such as biological systems that can change properties over time [19]. The NMJ modeled here is a particularly unique example of a dynamic biological system. The model includes a single release of acetylcholine (2000 molecules) from the neuron into a 50 nm2 region of the junction, containing 25 acetylcholinesterase molecules (biologically, these are tetramers treated individually) on each side of the junction and 50 nicotinic acetylcholine receptors (nAChR). When an individual acetylcholine molecule interacts with either the enzyme or the receptor, the agents both change.

Results and discussion

The model described here permits small molecular agents (i.e. ACh, inhibitor, and activator) to travel through the neuromuscular junction and interact with proteins (i.e. AChE or nAChR), binding and dissociating according to their state. Each agent is a biological entity and the interactions between protein and small molecule are based upon experimentally determined rate constants. As with the real-world, this model is limited in that interactions can only occur when two criteria are satisfied: two agents must be in physical proximity to each other, and they each must be capable of binding (i.e. no partner for ACh, and 0 or 1 partner for nAChR). To maintain consistency with the reality of a single endosome release, the model is spatially constrained along the y-axis to the junction distance, while the x-axis is allowed to remain unconstrained to simulate diffusion into and out of the region of interest.

To measure model performance and ensure the accuracy of the predictions generated, a set of simulations were performed wherein neither inhibitors nor activators were present and also with full AChE inhibition (Figure 2). The predicted shape and response to AChE inhibition is in agreement with in vitro observations [20, 21], and the predicted end-plate potential (e.p.p.; represented by number of open receptor channels in our model) rise-time of 35 msec and half-life of 94 msec are comparable with observed values of 6 msec and 7 msec, respectively [20]. The nearly 10-fold discrepancy between prediction and observation is not unexpected, as these references provide data from an amphibian leg muscle measurement and also demonstrated that both the e.p.p. rise-time and half-life can be affected by extracellular ionic strength, membrane potential, amount of ACh released, and the placement of electrodes for measurements. Previous ex vivo data in frog muscles show similar time scales and half-lives in the presence and absence of AChE inhibitors [21].Inhibition of AChE activity sharply decreases ACh clearance and elongates the e.p.p. (Figure 2). When maximally inhibited by GB (GB to AChE ratio of 2), the e.p.p. rise time was still 35 msec, however the half-life increased to 5.8 s. Inhibition at an inhibitor-enzyme ratio of 1 resulted in the appearance of two ACh turnover phases, the fast phase corresponding to the first 220 msec with a plateau equivalent to the fully inhibited level and a half-life of 2.3 msec, and a second, slower phase with a half-life of 211 msec. The e.p.p. prediction (Figure 2B) in this case also exhibited a two state response, with a signal plateau of 209 msec before signal decay with a half-life of 206 msec. Clearly, AChE inhibition affects multiple aspects of ACh clearance and the resulting e.p.p. duration.

Figure 2
figure 2

Model validation with active and GB inhibited AChE at two GB:AChE ratios (1, 2). Rates of ACh turnover are consistent with published values (A) and inhibition elongates the e.p.p. (B). Legend refers to both A and B.

To test the effects of general activation as a therapy for AChE inhibition by OPs, the impacts of small changes in activation parameters and activator concentrations on the e.p.p. were evaluated (Figure 3). By altering these parameters by 10, the inhibitor binding coefficient (ε-inh), was identified as the most critical parameter in affecting e.p.p. duration. Following the demonstration that modulating inhibitor binding could alleviate OP inhibition, the model was used to predict optimal activation coefficients to rescue severely inhibited NMJs. Considering this worst-case situation in which the majority of AChE would be inhibited, the inhibitor to activator ratio (IAR) at which the e.p.p. duration was 25% that of the uninhibited duration (Table 1) for several OP inhibitors was found, and used with a fully inhibited enzyme (inhibitor to enzyme ratio of 2:1) to calculate the lowest ε-inh for again returning the e.p.p. duration to within 25% of the uninhibited value. We hypothesize that the activator used here is allosterically activating, in that we observe changes in OP binding as being more critical than altering active site chemistry. Allosteric binding sites can be associated with long-range structural rearrangements [22] which could modulate the AChE active site cleft and “neck” without significantly altering the arrangement of the catalytic triad.

Figure 3
figure 3

Prediction of effects of activating AChE on e.p.p. duration in the presence of a 1:1 GB:AChE ratio. Varying activation parameters has a greater effect than varying activator concentration.

Table 1 Predicted therapeutic parameters for preventing OP-induced intoxication

By using an ABM, we take into account the dynamic nature of the activation event. In this case, the binding of a single activator agent to a single enzyme agent will transiently alter the behavior of that enzyme (for example, reducing the inhibitor’s binding affinity). With an ODE-based model, this small nuance of a single agent changing activity would be lost due to the averaging achieved by considering the populations of all enzyme agents. Within the context of this model, and the example just provided, a single activator binding to a single enzyme would account for a change of only 2% of the enzymes present, a negligible signal in the ODE yet one that we have shown to produce significant emergent behavior using the ABM.


Our model supports the hypothesis that as a potential therapeutic route, allosteric AChE activators provide a novel and useful approach for treating OP intoxication. Although our model does not rule out other methods and mechanisms of action, we contend that isolation and targeting the most effective allosteric therapies could produce life-saving effects by partially ameliorating the deleterious actions of AChE inhibitor binding.


Model development began with conceptually representative reaction diagrams, considering the binding of acetylcholine (ACh) to the nicotinic acetylcholine receptors (nAChR), ACh turnover by AChE, inhibition of AChE, and activation of AChE. These events can be represented by Figure 4, where A is ACh; R is the receptor; AR and AO are monoliganded closed and open receptors, respectively; A2R and A2O are diliganded closed and open receptors, respectively; A2D is the desensitized receptor; E is AChE; I is the inhibitor; and ε are the activation coefficients for substrate binding (ε-on and ε-off), catalysis (ε-cat), and inhibitor binding (ε-inh).

Figure 4
figure 4

Reaction diagram for nAChR-mediated signaling at the NMJ.

The system was modeled in NetLogo [23] as a single quantal injection of substrate [24]. The rules assigned to the individual agents are summarized below:

nAChR: cannot move; can bind one or two ACh molecules; can convert from closed and bound to open and bound with either one or two ACh ligands; can convert from open with two ligands to desensitized with two ligands; can release either ligand

AChE: cannot move; can bind single ACh molecule; can convert ACh to product; can bind inhibitor; can be aged by inhibitor; can bind activator; can release activator or ACh

ACh: can move with diffusion; can bind to nAChR or AChE; can be converted to product

Inhibitor: can move with diffusion; can bind to AChE; can induce AChE aging

Activator: can move with diffusion; can bind to AChE; can alter AChE activities

Product: disappears upon formation; increases accounting tracker by 1 unit

All of the microscopic rate constants were obtained directly from experimental reports, and there was no fitting of the model to observed data. The behavior resulting from each simulation run therefore emerges from the biophysical properties of the system and the definition of the rules. Receptor binding rates were parameterized from Hatton [25]. Enzyme turnover rates were obtained from Salih [26] and Hasinoff [27]. Inhibitor parameters were obtained from Hodge [28], Radic [29] and Worek [30] (Table 2).

Table 2 Inhibitor parameters for model

Supporting information available

The complete NetLogo code is included as supplemental information and is also available at the NetLogo model library ( Additional file 1.

Authors’ information

RRC is a protein biochemist, and an expert in enzyme inhibition. PJR is a biological modeler with extensive experience in ODE modeling in advanced biology, physics and mathematics. JJS is a board-certified toxicologist, a medicinal chemist, and expert in acetylcholinesterase inhibition by OPs. JMG is a board-certified toxicologist, a biological modeler, and expert in the toxicology and modeling of OP pesticides and nerve agents.


  1. Singh SP, Aggarwal AD, Oberoi SS, Aggarwal KK, Thind AS, Bhullar DS, Walia DS, Chahal PS: Study of poisoning trends in north India – a perspective in relation to world statistics. J Forensic Leg Med. 2013, 20 (1): 14-18. 10.1016/j.jflm.2012.04.034.

    Article  CAS  PubMed  Google Scholar 

  2. Jett DA, Yeung DT: The CounterACT Research Network: basic mechanisms and practical applications. Proc Am Thorac Soc. 2010, 7 (4): 254-256. 10.1513/pats.201001-003SM.

    Article  PubMed Central  PubMed  Google Scholar 

  3. Gupta B, Sharma R, Singh N, Kuca K, Acharya JR, Ghosh KK: In vitro reactivation kinetics of paraoxon- and DFP-inhibited electric eel AChE using mono- and bis-pyridinium oximes. Arch Toxicol. 2014, 88 (2): 381-390. 10.1007/s00204-013-1136-z.

    Article  CAS  PubMed  Google Scholar 

  4. Rosenberg YJ, Gearhart J, Mao L, Jiang X, Hernandez-Abanto S: Protection against paraoxon toxicity by an intravenous pretreatment with polyethylene-glycol-conjugated recombinant butyrylcholinesterase in macaques. Chem Biol Interact. 2014, 210: 20-25.

    Article  CAS  PubMed  Google Scholar 

  5. Nachon F, Brazzolotto X, Trovaslet M, Masson P: Progress in the development of enzyme-based nerve agent bioscavengers. Chem Biol Interact. 2013, 206 (3): 536-544. 10.1016/j.cbi.2013.06.012.

    Article  CAS  PubMed  Google Scholar 

  6. Abbas R, Hayton WL: A physiologically based pharmacokinetic and pharmacodynamic model for paraoxon in rainbow trout. Toxicol Appl Pharmacol. 1997, 145 (1): 192-201. 10.1006/taap.1997.8168.

    Article  CAS  PubMed  Google Scholar 

  7. Gearhart JM, Jepson GW, Clewell HJ, Andersen ME, Connolly RB: Physiologically based pharmacokinetic model for the inhibition of acetylcholinesterase by organophosphate esters. Environ Health Perspect. 1994, 102 (suppl. 11): 51-60.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  8. Poet TS, Kousba AA, Dennison SL, Timchalk C: Physiologically based pharmacokinetic/pharmacodynamic model for the organophosphorous pesticide diazinon. Neurotoxicology. 2004, 25 (6): 1013-1030. 10.1016/j.neuro.2004.03.002.

    Article  CAS  PubMed  Google Scholar 

  9. Timchalk C, Poet TS: Development of a physiologically based pharmacokinetic and pharmacodynamic model to determine dosimetry and cholinesterase inhibition for a binary mixture of chlorpyrifos and diazinon in the rat. Neurotoxicology. 2008, 29 (3): 428-443. 10.1016/j.neuro.2008.02.004.

    Article  CAS  PubMed  Google Scholar 

  10. Gearhart JM, Jepson GW, Clewell HJ, Andersen ME, Conolly RB: Physiologically based pharmacokinetic and pharmacodynamic model for the inhibition of acetylcholinesterase by diisopropylfluorophosphate. Toxicol Appl Pharmacol. 1990, 106 (2): 295-310. 10.1016/0041-008X(90)90249-T.

    Article  CAS  PubMed  Google Scholar 

  11. Chen K, Seng KY: Calibration and validation of a physiologically based model for soman intoxication in the rat, marmoset, guinea pig and pig. J Appl Toxicol. 2012, 32 (9): 673-686. 10.1002/jat.1671.

    Article  CAS  PubMed  Google Scholar 

  12. Van der Merwe D, Brooks JD, Gehring R, Baynes RE, Monteiro-Riviere NA, Rivierre JE: A physiologically based pharmacokinetic model of organophosphate dermal absorption. Toxicol Sci. 2006, 89 (1): 188-204.

    Article  CAS  PubMed  Google Scholar 

  13. Bonabeau E: Agent-based modeling: methods and techniques for simulating human systems. Proc Natl Acad Sci USA. 2002, 99 (suppl. 3): 7280-7287.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  14. An G, Mi Q, Dutta-Moscato J, Vodovotz Y: Agent-based models in translational systems biology. Wiley Interdiscip Rev Syst Biol Med. 2009, 1 (2): 159-171. 10.1002/wsbm.45.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  15. Pogson M, Holcomb M, Smallwood R, Qwarnstrom E: Introducing spatial information into predictive NF-kB modelling – an agent-based approach. PLoS ONE. 2008, 3 (6): e2367-10.1371/journal.pone.0002367.

    Article  PubMed Central  PubMed  Google Scholar 

  16. An G: A model of TLR4 signaling and tolerance using a qualitative, particle-event-based method: introduction of spatially configured stochastic reaction chambers (SCSRC). Math Biosci. 2009, 217 (1): 43-52. 10.1016/j.mbs.2008.10.001.

    Article  CAS  PubMed  Google Scholar 

  17. Zhang L, Athale CD, Deisboeck TS: Development of a three-dimensional multiscale agent-based tumor model: simulating gene-protein interaction profiles, cell phenotypes and multicellular patterns in brain cancer. J Theor Biol. 2007, 244: 96-107. 10.1016/j.jtbi.2006.06.034.

    Article  CAS  PubMed  Google Scholar 

  18. Kaul H, Ventikos Y: Investigating biocomplexity through the agent-based paradigm. Brief Bioinform. 2013, E-pub doi:10.1093/bib/bbt077

    Google Scholar 

  19. Lipniacki T, Puszynski K, Paszek P, Brasier AR, Kimmel M: Single TNFalpha trimmers mediating NF-kappaB activation: stochastic robustness of NF-kappaB signaling. BMC Bioinform. 2007, 8: 376-10.1186/1471-2105-8-376.

    Article  Google Scholar 

  20. Katz B, Miledi R: The binding of acetylcholine to receptors and its removal from the synaptic cleft. J Physiol. 1973, 231 (3): 549-574.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  21. Kordas M: On the role of junctional cholinesterase in determining the time course of the end-plate current. J Physiol. 1977, 270: 133-150.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  22. Sagermann M, Chapleau RR, DeLorimier E, Lei M: Using affinity chromatography to engineer and characterize pH-dependent protein switches. Prot Sci. 2009, 18 (1): 217-228.

    CAS  Google Scholar 

  23. Wilenski U: NetLogo. 1997, Evanston, IL: Center for Connected Learning and Computer-Based Modeling. Northwestern University, Available:

    Google Scholar 

  24. Land BR, Salpeter EE, Salpeter MM: Acetylcholine receptor site density affects the rising phase of miniature endplate currents. Proc Natl Acad Sci USA. 1980, 77 (6): 3736-3740. 10.1073/pnas.77.6.3736.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  25. Hatton CJ, Shelley C, Brydson M, Beeson D, Colquhoun D: Properties of the human muscle nicotinic receptor, and of the slow-channel myasthenic syndrome mutant εL221F, inferred from maximum likelihood fits. J Physiol. 2003, 547 (3): 729-760. 10.1113/jphysiol.2002.034173.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  26. Salih E: Catalysis by acetylcholinesterase in two-hydronic-reactive states: integrity of deuterium oxide effects and hydron inventories. Biochem J. 1992, 285: 451-460.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  27. Hasinoff BB: Kinetics of acetylthiocholine binding to electric eel acetylcholinesterase in glycerol/water solvents of increased viscosity: Evidence for a diffusion-controlled reaction. Biochim Biophys Acta. 1982, 704 (1): 52-58. 10.1016/0167-4838(82)90131-5.

    Article  CAS  PubMed  Google Scholar 

  28. Hodge AS, Humphrey DR, Rosenberry TL: Ambenonium is a rapidly reversible noncovalent inhibitor of acetylcholinesterase, with one of the highest known affinities. Mol Pharmacol. 1992, 41 (5): 937-942.

    CAS  PubMed  Google Scholar 

  29. Radic Z, Taylor P: Interaction kinetics of reversible inhibitors and substrates with acetylcholinesterase and its fasciculin 2 complex. J Biol Chem. 2001, 276 (7): 4622-4633. 10.1074/jbc.M006855200.

    Article  CAS  PubMed  Google Scholar 

  30. Worek F, Thiermann H, Szinicz L, Eyer P: Kinetic analysis of interactions between human acetylcholinesterase, structurally different organophosphorus compounds and oximes. Biochem Pharmacol. 2004, 68 (11): 2237-2248. 10.1016/j.bcp.2004.07.038.

    Article  CAS  PubMed  Google Scholar 

Download references


The authors thank Dr. Yaroslav Chushak, C. Erik Hack, Tammie Covington, Christopher Ruark and Amanda Hanes for comments and insight.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Richard R Chapleau.

Additional information

Competing interests

The authors declare that they have no competing financial interests.

Authors’ contributions

RRC designed the scheme, constructed the model, performed the simulations, analyzed the data, and drafted the manuscript. PJR designed the scheme, reviewed the model, and drafted the manuscript. JJS reviewed the data and drafted the manuscript. JMG parameterized and reviewed the model, directed the simulations, analyzed the data and drafted the manuscript. All authors read and approved the final manuscript.

Electronic supplementary material


Additional file 1: Supporting Information Legends. NetLogo Code. Text code for the model. This code functions, but does not include the graphical user interface for the model, which is available at the NetLogo model library. (DOCX 35 KB)

Authors’ original submitted files for images

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Chapleau, R.R., Robinson, P.J., Schlager, J.J. et al. Potential new therapeutic modality revealed through agent-based modeling of the neuromuscular junction and acetylcholinesterase inhibition. Theor Biol Med Model 11, 42 (2014).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Neuromuscular junction
  • Acetylcholinesterase
  • Organophosphate intoxication
  • Medical countermeasure
  • Allosteric enzyme regulation
  • Agent-based model