Eukaryotic protein kinases (PKs) are a large family of proteins critical for cellular response to external signals, acting as molecular switches. PKs propagate biochemical signals by catalyzing phosphorylation of other proteins, including other PKs, which can undergo conformational changes upon phosphorylation and catalyze further phosphorylations. Although PKs have been studied thoroughly across the domains of life, the structures of these proteins are sparsely understood in numerous groups of organisms, including plants. In addition to efforts towards determining crystal structures of PKs, research on human PKs has incorporated molecular dynamics (MD) simulations to study the conformational dynamics underlying the switching of PK function. This approach of experimental structural biology coupled with computational biophysics has led to improved understanding of how PKs become catalytically active and why mutations cause pathological PK behavior, at spatial and temporal resolutions inaccessible to current experimental methods alone. In this review, we argue for the value of applying MD simulation to plant PKs. We review the basics of MD simulation methodology, the successes achieved through MD simulation in animal PKs, and current work on plant PKs using MD simulation. We conclude with a discussion of the future of MD simulations and plant PKs, arguing for the importance of molecular simulation in the future of plant PK research.
Eukaryotic cells are able to rapidly respond to biotic and abiotic changes in their environment, allowing them to follow chemical gradients , alter their growth cycle in response to changes in chemical signal concentrations , or co-ordinate their behavior with other cells , among numerous other abilities. Though the mechanisms for transduction of external signals into cellular response vary in detail across different organisms, cell types, and signal types, eukaryotic protein kinases (PKs) are central to many of the responsible signaling pathways . The proximate function of PKs is to catalyze the transfer of phosphoryl groups, typically donated by ATP molecules, to the Tyr, Thr, or Ser residues of other proteins. Phosphorylation can alter the physical and catalytic properties of proteins, including PKs themselves. The propagation of a signal through successive phosphorylation of PKs by other PKs is often referred to as a kinase cascade and is an important component of the eukaryotic cellular circuitry . Ultimately, kinase cascades are often coupled to gene regulation machinery, completing the connection between the environment of a cell and its gene expression.
As a consequence of their importance in such a wide variety of cellular processes, PKs with abnormal function or expression levels are often implicated in disease, perhaps most notably in many cancer types [6–8]. For this reason, PKs attract considerable attention from researchers in a biomedical context. A large portion of this effort has been directed toward understanding how PKs are structured and how their structures change upon phosphorylation or interaction with other proteins . Owing to their primary role in information relay, PKs have evolved to act as highly sensitive switches rather than catalytic workhorses rapidly churning out products , meaning that understanding the structure and conformational dynamics of PKs is central to understanding their function and dysfunction within cells. Further supporting the importance of a structural understanding of PKs to medicine is the rise of rational PK inhibitor design using computational and experimental models, which depend on the detailed knowledge of protein structure [11,12].
In this review, we aim to highlight the importance of MD simulations for understanding the function and regulation of plant PKs. In the first section, we discuss the state of plant PK structural biophysics and describe general PK structure. Next, we give a brief overview of molecular dynamics (MD) simulation, emphasizing the scientific questions motivating its development. We describe past work on non-plant eukaryotic PKs and the unique biological insights gained through combining MD simulation and experimental results. We then review the few examples of MD simulation applied to plant PKs. Finally, we discuss the possible role of MD simulation in the study of plant PKs and argue that MD simulation, complementing structural experimental studies, will be critical to further progress in the field.
Plant PK structure and dynamics: where few have gone before
The basic anatomy of a typical PK is shown in Figure 1a, with some of the most important features highlighted. The αC helix, shown in gray, must be folded and positioned close to the remainder of the kinase domain for catalysis. The Lys-Glu bond is an ionic interaction between a conserved Lys-Glu pair, where the Glu is part of the αC helix, which must be formed in catalytically competent conformations . While ATP is not present in this structure, ATP must be bound at the labeled ATP-binding site with the Gly-rich loop closing over it, while the Phe in the DFG motif must be positioned outward, as in the displayed structure. Proper positioning of the Phe in the DFG motif and of the αC helix completes a contiguous stack of hydrophobic residues, known as the regulatory spine (R-spine) . The activation loop (A-loop) must be unfolded flat over the top of the C-lobe, allowing for substrate binding and exposing the catalytic residues in the catalytic loop (C-loop) to the substrate. Positioning the A-loop for activation is often accomplished through phosphorylation and accounts for the activation of PKs by other PKs.
Structural overview of plant protein kinases.
In contrast with human and, more generally, mammalian PKs, plant PK research has focused far less on structure and conformational dynamics; perhaps, the moral and financial imperatives driving human PK research at such a high level of spatio-temporal detail are not as readily evident for plants. Regardless of the reasons for these differences, plant PK research has made significant strides toward understanding PK function within a variety of signaling pathways exemplified by research on brassinosteroid signaling [20–24]. Still, the dearth of plant PK structures determined by any experimental method (11 entries in total), quantified by the number of unique Arabidopsis thaliana kinase domains with at least one structure on the Protein Data Bank (PDB) , has limited a detailed understanding of plant signaling pathways (Figure 1b). We analyzed the availability of PK structures from the model organism A. thaliana rather than any plant of any direct relevance to agriculture due to the even further scarcity of PK structures from any other plant, where we found two entries from Zea mays (maize) and one entry from Solanum pimpinellifolium (tomato).1
The present and projected disastrous effects of global climate change on agricultural growing conditions in regions with high levels of food insecurity , coupled with a growing global population predicted to outgrow current food production projections , make the need for crop innovation increasingly clear. The importance of PKs in plant growth and drought resistance could make them prime targets for genetic engineering or agrochemical development . For example, changes in the expression of several PKs have been found to improve plant drought resistance [28–31], while expression of constitutively active PKs in plants has been shown to improve stress tolerance , increase sensitivity to stress signaling hormones , and improve leaf growth . A detailed understanding of PK structure and dynamics is invaluable for guiding the introduction of point mutations or more extensive modifications to PK sequences in order to tune function and ultimately improve crop yield. Just as structural information guides rational design of inhibitors for mammalian PKs involved in cancers [11,12], agrochemicals could be rationally designed to target plant PKs given sufficient knowledge of their structure, providing another avenue for the manipulation of plant growth and stress resistance.
The study of plant PKs from a structural–dynamics perspective also has the potential to uncover previously unknown modes of PK regulation. Considering the extensive growth of plant kinomes through the history of evolution partially owing to adaption to environmental stresses [35,36] in addition to the broad repertoire of post-translational modifications (PTMs) utilized by plants to alter protein function [37,38], it seems likely that unique modes of PK regulation await discovery within the kinomes of plants. Additionally, a structural understanding of how PKs fit into signaling pathways, in terms of their protein–protein interactions, is critical for insights into the biological function of PKs in the cellular context. For example, a dimeric crystal structure was found for subfamily 2 Snfl-related kinase 6 (SnRK2.6) and a type 2C protein phosphatase (HAB1) in A. thaliana . The binding interface of HAB1 and SnRK2.6 closely resembles that of HAB1 and the abscisic acid (ABA) receptor PYL2 when ABA is bound, illuminating how ABA signaling is coupled with SnRK2.6 enzymatic control. With structural models of protein–protein interactions involving PKs, it could be possible to rationally engineer PKs through mutagenesis of interfacial residues to specifically target proteins in one downstream signaling pathway over another, providing fine control over cellular signaling. One potential target for this mode of study is BRI1-ASSOCIATED KINASE 1 (BAK1), a receptor-like kinase involved in both plant growth and immune signaling, where it is already known that differential phosphorylation of BAK1 leads to differential specificity in its phosphotransfer activity . As is the case with many plant PKs, the structural details needed to explain the underlying mechanisms of these changes in specificity are missing.
If we are to understand the structural and dynamic basis for plant PK function, how can the successes of the human PK research program be translated into advances in plant biology? While the human kinome is far from fully characterized on a structural–dynamics level, great progress has been made by combining experimental methods, such as X-ray crystallography and nuclear magnetic resonance (NMR) spectroscopy, with computational methods such as MD simulation. While the information provided by X-ray crystallography primarily concerns the static ‘native’ structure of a protein, NMR experiments can provide the timescales of interconversion between different protein conformations along with some structural details. Additionally, computational structure prediction methods, ranging from homology modeling using the extensive library of protein folds in experimentally solved structures [41–44] to physics-based modeling guided by evolutionary information [45–49], have great potential to extend static structural knowledge of plant kinomes and compensate for the lack of experimental structures. MD simulation, however, enables researchers to gain unique insights into the movement of each atom in a protein over relatively long periods of time (nanoseconds to milliseconds).
Molecular dynamics simulation
The basics of classical MD simulations
MD simulation was initially developed by Fermi et al.  in order to study the dynamics of relatively simple multiparticle systems with no closed-form mathematical solutions. From then on, MD simulations of gaseous argon , liquid water , and eventually proteins  were performed as the required software and hardware improved, and as protein crystal structures began to be solved . MD simulation is often likened to a ‘computational microscope’ [55,56] in that it can give a detailed, atom-by-atom view of how proteins and nucleic acids move with time .
Modern MD simulation of proteins requires several sources of input. First, an initial set of protein co-ordinates is needed from which to start simulations. These initial co-ordinates are typically taken from X-ray crystallography experiments or are inferred through homology modeling, which uses proteins with similar sequences and known structures to predict a possible structure. In principle, one could construct a computational model of an unfolded protein from sequence information alone and use MD simulation to observe the folding process [58,59], but this approach is usually impractical due to the enormous amount of computing time required to fold most proteins of biological interest.
The second input needed for MD simulation is the set of mathematical functions defining forces between atoms in the protein. One must decide both the forms of the functions used to calculate forces, approximating more accurate physical theory, as well as the parameters tuning those functions to match with experimental quantities and quantum calculations (Figure 2). In practice, for simulation of a protein consisting of the 20 standard amino acids, one must only choose from one of the many available biomolecular force fields [60–63], which provide a functional form for force calculations as well as a set of parameters carefully tuned to reproduce certain experimental results.
Interatomic interactions in classical MD force fields.
Once the initial atomic co-ordinates have been acquired and the force field has been chosen, MD simulation can begin.2 At their core, MD programs operate using the same physical principles that Isaac Newton developed over 300 years ago to describe the motion of far larger objects, hence the word ‘classical’ in the title of this section. At the same time, methodological advances continue to improve the accuracy, efficiency, and capability of MD simulation .
The motion of each atom is calculated in time steps, producing a sequence of protein conformations describing its motion with time. To give a simple illustrative example of the process, we consider a single atom moving along a straight line (Figure 3a). The position of the atom on the line at a certain time is denoted as x(t) with a velocity v(t) and acceleration a(t). Just as one might calculate the motion of an object in a secondary school physics problem, the position of the atom after a certain amount of time, which we call Δt, can be calculated according to3x(t + Δt) ≈ x(t) + v(t)Δt + (1/2m)f(t)Δt2. The initial position and velocity of the particle are provided at the start of the simulation. From Newton's second law of motion, f = ma, we know that the acceleration of an object is proportional to the force that it feels, meaning that the force fields described above will provide the a(t) term [where a(t) = f(t)/m]. However, our system has only one particle with no forces acting on it and will perpetually continue motion in one direction (Figure 3b), which is not particularly useful for our purposes. If we add another atom so that the two atoms are attached by a spring and the force between them is f(x1, x2) = −k|x1 − x2| (Figure 3c), where | · | is the absolute value of ‘·’, we will have to calculate the motion of both atoms, taking into account the force that they impart on one another (Figure 3d). If we instead consider a system with many interacting atoms each moving in three dimensions, the dynamics can become very complicated and the number of force calculations between each pair of atoms will increase rapidly with the total number of atoms in the system.
Two examples of simple systems with dynamics calculated using Newtonian mechanics.
Analysis and extensions of MD
Once trajectories are completed, useful information has to be extracted in order to provide biological insight. One common method of analyzing MD trajectories is to calculate a free energy landscape (FEL) over protein structural features (Figure 4). A FEL provides information about the relative stability of different protein conformations as well as the free energy barriers between stable states, which control how fast a protein can switch between different conformations. For example, if two conformations of a protein are compared, one can use the free energy difference between the two to determine their relative stabilities which can be translated into the relative amount of time a given protein will spend in either conformation. A glaring limitation of the FEL is that proteins require a large number of variables to fully describe the motion of all atoms therein; one must choose several ‘important’ variables to calculate an FEL which can be visualized. Deciding what makes a variable important and how to identify such variables is a challenge and remains an active area of research.
Timescales in protein dynamics.
Beyond the basic MD simulations we have described above, often referred to as unbiased MD for reasons which will become clear shortly, many other computational methods for studying protein conformational dynamics exist. One class of these methods is biased MD simulation, where the potential energy function determining the motion of a protein is altered in some way from the realistic approximate energy function of unbiased MD, promoting conformational changes that could otherwise take prohibitively long times to occur . Because the bias added to the potential energy function is a known quantity, one can account for the bias in the analysis stage and produce the same FEL that might have come from unbiased simulation, but with less computational time required. Despite the potential for increased efficiency in biased simulations, it can be a challenge to determine how exactly to add the biasing energy to a simulation in order to actually improve efficiency . Another method that can be used involves creating coarse-grained models of a protein structure, where multiple atoms are combined into representative units and the motion of these units is calculated [66–68]. This strategy provides several advantages, decreasing the number of calculations a computer must do and thus allowing simulation of larger systems for longer periods of time, while sacrificing the detailed description of dynamics provided by all-atom MD .
MD simulation and structural biology: a successful approach for animal PKs
We have discussed the general benefits of MD simulation in the previous section, but how has MD contributed to the study of PKs? For animal PKs, an approach combining structural biology experiments with MD simulation has uncovered the dynamic nature of PK structure for many different specific proteins . MD simulations of PKs have uncovered the effects of phosphorylation and ATP binding, including induction of conformational change along activation pathways [70–77], functional mechanisms of activating mutations [78–80], binding properties of inhibitors [81–85], as well as allosteric interactions between different structural domains [86,87]. All of these properties are difficult to access experimentally but are critical for fully understanding, and eventually manipulating, PKs as the molecular switches that they are. It is precisely because of these unique abilities that MD simulation has played such an important role in animal PK research, and will very likely continue to do so for years to come. Although MD simulation provides unique insights into PK dynamic behavior, the successes of MD as a method of studying PKs have come through integration with experimental results. Information can flow either way; MD simulation is capable of providing a wide range of experimentally testable predictions, while experimental results can be used to improve the accuracy and efficiency of simulations. For example, observables from NMR spectroscopy experiments can be compared with the same quantities calculated from MD simulations or used to restrain atoms in simulations in such a way as to ensure agreement with experimental values [88,89]. The idea behind introducing NMR restraints is to account for inaccuracies in the empirical force fields, with the intent of improving the predictive power of simulations for quantities not included in the restraints. As opposed to the quantitative physical properties measured in NMR, MD simulations can also be used to make qualitative predictions, such as which residues are important for protein function or stability [90–92]. This type of prediction can be tested through site-directed mutagenesis coupled with functional assays. Inversely, deep mutational scanning  can be used to experimentally identify functionally important residues in a protein, and MD simulation can be used to understand why a particular residue is important .
Activation pathways and post-translational modifications
Out of the PK properties typically studied by MD simulation, the activation pathway, the most probable sequence of conformations taken to travel from the inactive to active state, is perhaps of the most fundamental importance. Characterizing how phosphorylation or release of an inhibiting domain leads to PK activation is analogous to understanding how gears and levers allow the flip of a mechanical switch to perform its ultimate function. However, PK activation pathways are difficult to characterize due to the long timescales over which conformational changes take place, meaning that large amounts of computer time are needed to observe an activation event, let alone other modes of PK motion [71,95].
Still, many groups have had success either through clever use of enhanced sampling techniques or through brute force simulation. Shukla and co-workers  used extensive MD simulations performed on the Folding@home distributed computing platform  to study the activation pathway of human c-Src kinase, describing the conformational path taken between crystal structures of the active and inactive states (Figure 5). While this transition is driven in vivo by phosphorylation, the authors utilized a sophisticated algorithm, called the string method in collective variables [97,98], in order to identify a likely transition path for unphosphorylated c-Src. The structures along this pathway were then used to start simulations which were, in turn, used to create a Markov state model (MSM). MSMs are approximate kinetic models created by grouping structurally similar conformations into discrete states and calculating rates of interconversion between states from the original MD trajectories [99,100]. Through this approach, they were able to predict the presence of intermediate conformational states, one of which was found to be similar to an inhibitor-bound X-ray crystal structure of cyclin-dependent kinase 2 (CDK2) . Further simulations of c-Src bound to the CDK2 inhibitor revealed that the inhibitor stabilized the c-Src inactive state, demonstrating the presence of an inhibitor-binding pocket present in the conformations between the active and inactive states.
Activation of human c-Src kinase.
The study of PK response to PTMs and protein–protein interaction is intimately linked to the study of activation pathways and is another integral part of the structural research on PKs. Past research has focused on the effects of PTMs or interacting proteins on both the protein conformational FEL [74–76] and directly on the efficiency of the phosphoryl transfer reactions catalyzed by PKs . In an example of the former case, Kuzmanic and co-workers examined the effects of A-loop phosphorylation and ATP binding to the human PK p38α using an enhanced MD sampling method called parallel tempering metadynamics [76,102]. p38α was simulated when unphosphorylated, doubly phosphorylated, doubly phosphorylated with ATP bound, and doubly phosphorylated with ATP and a docking peptide bound in order to study the effects of these factors on the p38α FEL. They found that in order for a stable state representing the active-like crystal structure to occur, phosphorylation, ATP binding, and docking peptide binding must all occur, even though only phosphorylation is present in the crystal structure. However, this apparent contradiction was resolved by showing that in simulations of a crystal super cell, the crystal contacts stabilized the active-like state leading to the discrepancies between X-ray crystal structures and MD simulation together with prior NMR results.
Given the high levels of PK structural conservation across the eukaryotes , we expect that the insights into PK activation gained from MD simulation of animal PKs to be of great use from understanding plant PK activation. The principles of R-spine alignment, including αC helix placement and DFG positioning, and A-loop unfolding together with other conserved motifs mentioned above, can be used as starting points for the study of plant PKs, although the possibility for new regulatory motifs remains enticing given the large size of plant kinomes .
Understanding PK dysfunction: oncogenic mutations and drug design
While a mechanistic understanding of normal PK function is of great scientific value, the importance of abnormally functioning PKs in cancer brings the focus of biomedical research toward understanding how mutations lead to pathological PK function and how inhibitors might be designed to treat specific cancers by blocking the function of overactive PKs. Through MD simulations, researchers can correlate differences in residue interactions upon mutation with shifts in conformational dynamics in order to understand how mutations alter PK function. Shan et al.  took this approach to study how an oncogenic mutation in the human epidermal growth factor receptor (EGFR) kinase domain leads to pathogenic activation. EGFR is a transmembrane receptor tyrosine kinase which forms homodimers upon binding of EGF , although the L834R mutation in the kinase domain causes constitutive EGFR activation . Using MD simulations validated by biochemical experiments, Shan and co-workers showed that the αC helix of the wild-type EGFR kinase domain displayed a high degree of disorder in monomeric form, while dimerization suppressed this disorder leading to activation. Simulation of the L834R mutant revealed suppressed disorder in the αC helix, promoting dimerization and therefore activation of EGFR.
Oncogenic PKs are a common target for the development of cancer drugs . Numerous studies have used MD simulation in order to characterize inhibitor binding to PKs, including attempts to explain differences in inhibitor affinity to different PKs [81,83,106] and to discover binding pathways . Binding pathways, in particular, are uniquely accessible by MD simulation, where alternative binding sites and important protein residue–inhibitor interactions can be observed . Both approaches to understanding the impact of cancer-causing mutations on animal PKs and designing inhibitory drugs could be extended to genetic engineering of plant PKs and the development of agrochemicals targeting plant PKs. In terms of agrochemical development, the established field of plant chemical genetics , where the focus of research is on screening small molecules for biological activity and identifying target proteins, could be complemented by the knowledge of plant PK structure and dynamics. Aside from fundamental differences in the goals of applying methodologies previously used to understand and treat cancer (increasing rather than inhibiting growth), the basic principles appear almost wholly transferrable to problems in plant cellular signaling and crop engineering.
MD simulation and structural biology: a successful approach for plant PKs?
Despite the successes achieved for animal PKs through marrying structural experiments and MD simulation, there has been relatively little progress on either front for plant PKs. Of the 940 known PKs in the A. thaliana kinome, only 27 have an entry in the PDB, 11 of which contain any part of a kinase domain.4 Compared with the human kinome, the structural features of plant kinomes are largely uncharacterized (Figure 1b), while details of conformational dynamics are practically unknown. In this section, we review what progress has been made on understanding the structure and dynamics of plant PKs.
In terms of molecular simulation, very little work was focused on plant PKs, let alone plant proteins in general. Two studies [108,109] have focused on calculating binding enthalpy of members of the brassinosteroid plant hormone group to their primary receptors, BRASSINOSTEROID-INSENSITIVE 1 (BRI1) and BAK1, both leucine-rich repeat receptor-like kinases. As brassinosteroids bind to the extracellular domains of BRI1 and BAK1 across the cell membrane from the cytoplasmic kinase domains (Figure 6), these studies provide important insights into brassinosteroid signaling but are of little relevance to the study of the plant PKs. As far as we are aware, the only MD simulation studies, to date, on plant PK kinase domains are our own, focusing on the conformational dynamics of BRI1 and BAK1 [110,111]. In our first paper, we sought to characterize the long-term conformational dynamics of the BRI1 and BAK1 kinase domains in their active phosphorylation states. In a second paper, we examined the structural mechanism of BAK1 deactivation by S-glutathionylation, which was observed from in vitro experiments . The remainder of this section will be dedicated to reviewing the results of these two papers.
A cartoon representation of the current state of structural knowledge for brassinosteroid signaling initiation.
BRI1 and BAK1 conformational dynamics
Despite the importance of brassinosteroid signaling through BRI1 and BAK1 for plant growth and development [40,113,114], the atomic-level details of early BRI1/BAK1 signaling events are largely unknown. In the generally accepted model of brassinosteroid signaling initiation, BRI1 binds to a single BL molecule, promoting association with BAK1 and reciprocal phosphorylation of the BRI1 and BAK1 kinase domains (Figure 6) . In addition to kinase domain phosphorylation, the largely unstructured BRI1 kinase inhibitor 1 is phosphorylated and subsequently released . The activated BRI1 kinase domain then phosphorylates downstream targets, ultimately leading to activation of the transcription factors BRASSINOZOLE-RESISTANT 1 and BRI1-EMS-SUPPRESSOR 1  and alteration of gene expression controlling plant growth and development. While this qualitative model has been pieced together through years of experiments [20,117–119], many questions remain concerning the precise mechanisms of BRI1 and BAK1 activation. For example, it is currently unknown how the BRI1–BAK1 transmembrane helices and kinase domains interact with one another and how BL binding activates BRI1 and BAK1. Improved understanding of the molecular basis for plant growth signaling will not only mark a significant advance in fundamental plant science, but also could inform rational genetic engineering efforts and agrochemical development targeting brassinosteroid signaling pathways to improve crop yield and stress resistance.
To address this knowledge gap, we sought to characterize the conformational dynamics of the BRI1 and BAK1 kinase domains in ATP-bound and fully phosphorylated forms by performing MD simulations on both kinase domains separately. We found that much like p38α , phosphorylation and ATP binding were insufficient for either protein to maintain active-like structures similar to their crystal structures. For BAK1, we found that the αC helix had a high propensity for disorder (Figure 7a), much like EGFR , while the BRI1 αC helix adopted both disordered and Src-like  inactive conformations (Figure 7b). To test whether the high degree of BAK1 αC helix disorder was due to an intrinsically disordered sequence or the surrounding BAK1 residues, we substituted the stable αC helix of human c-Src into our BAK1 computational model and performed further MD simulations. The BAK1–Src chimera exhibited an enhanced level of αC helix stability, suggesting that intrinsic disorder in the native BAK1 αC helix sequence is at least partly responsible for its unfolding in simulations. We also performed circular dichroism experiments on the BRI1 αC helix peptide in order to characterize its propensity for disorder in isolation, finding that the BRI1 αC helix is less disordered than the EGFR αC helix, but more disordered than the c-Src αC helix.
BRI1 and BAK1 kinase domain dynamics.
Finding disordered αC helices in BRI1 and BAK1 led us to ask how common this feature is across entire plant kinomes. Accordingly, we used the PONDR-FIT disorder-prediction webserver  to predict disorder in putative αC helix regions of all 940 known A. thaliana PKs , finding that αC helix disorder could be a common attribute in plant kinomes (Figure 7c,d). As phosphorylation was insufficient in our simulations to keep the kinase overwhelmingly in an active-like state, we proposed that some other factor not included in our simulations may play a role in BRI1 and BAK1 activation by suppressing αC helix disorder, possibly the actual physical interaction between the two kinase domains. In any case, as αC helix disorder is predicted to be widespread in the A. thaliana kinome, it may be a widely used regulatory mechanism in plant PKs.
BAK1 deactivation by S-glutathionylation
In addition to activation by phosphorylation, BAK1 can be deactivated through S-glutathionylation . S-glutathionylation is a post-translational modification where glutathione, a ubiquitous tripeptide shown to be critical for normal function of several plant and animal species [122,123], forms a disulfide bond with a Cys residue. While Bender and co-workers identified three possible S-glutathionylation sites on BAK1, it was unclear which site or which combinations of sites were responsible for BAK1 deactivation as well as how, at the molecular level, S-glutathionylation was able to deactivate BAK1. To understand the mechanism of BAK1 deactivation by S-glutathionylation, we performed MD simulations of BAK1 singly modified at each possible site for a total of three simulation systems. Using the unglutathionylated BAK1 simulations as a reference, we investigated the effects of BAK1 S-glutathionylation on its FEL. We calculated the Kullback–Leibler divergence, a measure of similarity between probability distributions, between the distributions of residue conformations in S-glutathionylated and unglutathionylated BAK1 simulations, finding that modification of one site, C408, had a large effect on BAK1 conformation across the kinase domain (Figure 8a,b). Further supporting C408 as an inhibitory S-glutathionylation site is the shift of the BAK1 FEL away from the active-like structure with C408 S-glutathionylation, while modification of other sites yielded only slightly changed landscapes (Figure 8c). While C408 S-glutathionylation induced long-range allosteric effects across the BAK1 kinase domain, αC helix unfolding is promoted through direct interaction with glutathione.
Deactivation of BAK1 by S-glutathionylation.
The future of MD and the plant kinome
Toward the lofty goal of characterizing the dynamics of any large set of proteins, it will be necessary to remain on the cutting edge of molecular simulation methodology. While the analogy of a computational microscope paints a vivid picture of the capabilities of MD simulation, it is important not to take the analogy too literally since the accuracy of simulation results is limited by the error inherent to numerical simulation, the approximate nature of calculated interatomic forces, and the high computational cost necessary for simulations on biologically relevant timescales . Despite all of the advances in MD methodology, there is much room for improvement. Perhaps, most clear among possible strategies for improving MD simulation are further integration of experimental results into simulations [88,125], benefitting both the experimentalist and the computational scientist. Additionally, successes of coarse-grained  and quantum  molecular modeling applied to animal PKs, allowing the study of larger and smaller spatio-temporal scales than classical all-atom MD, respectively, demonstrate the potential of methodological developments in other simulation methods for the study of plant PKs.
For the use of physics-based computational models to find success in the study of plant PKs, and plant proteins in general, it will be necessary for future generations of plant biologists to be trained in quantitative thinking and methods . As we have pointed out, very little structural information exists concerning plant PKs despite their enormous number within plant species . Researchers with strong backgrounds in physics, mathematics, and computational science will be necessary to uncover the atom-level details of plant PKs and in doing so, guide future crop innovation. We emphasize that this directive should not be limited to the PKs; the study of plant proteins, in general, will greatly benefit from realization of the vision outlined by many leading plant biologists .
Flowering plants tend to have large kinomes, ranging from 600 to 2500 members among analyzed species . In contrast, the human kinome consists of ∼500 members . A possible explanation is that plants have such extensive kinomes due to their need to respond to constantly changing external stimuli without the benefits of motility or a central nervous system . Regardless of the reasons for the size of plant kinomes, it is clear that plant PKs are underrepresented in the PDB with 296 kinase domain entries from unique Homo sapiens proteins (55.0% of the human kinome) when compared with 27 entries from unique A. thaliana proteins (2.9% of the A. thaliana kinome), 2 entries from Z. mays, and 1 entry from S. pimpinellifolium. While the focus on human PKs is understandable due to their importance in disease, we argue that a research program is needed for plant PKs mirroring that of human PKs. First, on scientific grounds, plant kinomes represent fertile ground for new discoveries, given the large number of PKs with no structural or dynamic characterization. While the catalytically active conformations of PKs are necessarily similar to one another, the inactive conformations of animal PKs are far more diverse, reflecting the effectively infinite number of ways a protein can be catalytically incompetent for a particular reaction [121,128]. Second, on practical grounds, understanding how plant PKs function could greatly improve genetic engineering and rational agrochemical design attempts to improve crop growth or stress resistance. For both of these reasons, it is time for the era of plant PK structural and dynamic proteomics.
BRI1-ASSOCIATED KINASE 1
cyclin-dependent kinase 2
epidermal growth factor receptor
free energy landscape
Markov state model
nuclear magnetic resonance
Protein Data Bank
2 Snfl-related kinase 6
We thank Shukla research group members for their critical reading of this manuscript. We also thank Steven C. Huber for the discussion of the topics covered in this review. A.S.M. and D.S. acknowledge the University of Illinois at Urbana-Champaign for the support.
The Authors declare that there are no competing interests associated with the manuscript.
The number of unique entries in the PDB for plant PK domains is higher than reported here, as we noted that the enzyme classification for at least several entries was misannotated in downloaded reports, classifying protein kinases as non-protein kinases. For the sake of simplicity, we have reported the numbers according to these reports.
†For the sake of brevity, we avoid discussion of several crucial steps in between securing input information and actual MD simulation as they are not important for the overall description of MD simulation. At the same time, the omitted steps are oftentimes among the most challenging in the entire process and cannot be ignored in a more thorough discussion.
‡In practice, this simple algorithm for calculating trajectories is not used due to the high degree of error it incurs. This error arises due to the necessarily non-vanishing values of Δt used in computer simulations, where v(t) and a(t) are approximated as constants over the period of this time step.
§As mentioned earlier for maize and tomato PK structures, this number could be an underestimate if the misannotation of enzyme classifications is a recurring issue in the PDB.