In SilicoModel-driven Assessment of the Effects of Brain-derived Neurotrophic Factor Deficiency on Glutamate and Gamma-Aminobutyric Acid: Implications for Understanding Schizophrenia Pathophysiology.

OBJECTIVE: Deficient brain-derived neurotrophic factor (BDNF) is one of the important mechanisms underlying the neuroplasticity abnormalities in schizophrenia. Aberration in BDNF signaling pathways directly or circuitously influences neurotransmitters like glutamate and gamma-aminobutyric acid (GABA). For the first time, this study attempts to construct and simulate the BDNF-neurotransmitter network in order to assess the effects of BDNF deficiency on glutamate and GABA.
METHODS: Using CellDesigner, we modeled BDNF interactions with calcium influx via N-methyl-D-aspartate receptor (NMDAR)- Calmodulin activation; synthesis of GABA via cell cycle regulators protein kinase B, glycogen synthase kinase and β-catenin; transportation of glutamate and GABA. Steady state stability, perturbation time-course simulation and sensitivity analysis were performed in COPASI after assigning the kinetic functions, optimizing the unknown parameters using random search and genetic algorithm.
RESULTS: Study observations suggest that increased glutamate in hippocampus, similar to that seen in schizophrenia, could potentially be contributed by indirect pathway originated from BDNF. Deficient BDNF could suppress Glutamate decarboxylase 67-mediated GABA synthesis. Further, deficient BDNF corresponded to impaired transport via vesicular glutamate transporter, thereby further increasing the intracellular glutamate in GABAergic and glutamatergic cells. BDNF also altered calcium dependent neuroplasticity via NMDAR modulation. Sensitivity analysis showed that Calmodulin, cAMP response element-binding protein (CREB) and CREB regulated transcription coactivator-1 played significant role in this network.
CONCLUSION: The study presents in silicoquantitative model of biochemical network constituting the key signaling molecules implicated in schizophrenia pathogenesis. It provides mechanistic insights into putative contribution of deficient BNDF towards alterations in neurotransmitters and neuroplasticity that are consistent with current understanding of the disorder.

INTRODUCTION

Schizophrenia is considered to be one among the top ten leading cause of disease related disability1) with a developing risk of ~0.7%.2) The etiology and pathophysiology remains unknown despite of several years of study and research. Various models and hypothesis emerged to better understand this disorder.3) With regards to neurotransmitter models, dopamine aberrations contribute critically to the pathogenesis of this disorder—hyper-dopaminergia is postulated to underlie the genesis of positive symptoms4) and hypodopaminergia to underlie the genesis of negative symptoms.5) There dopamine abnormalities are considered to be downstream effects of cortical glutamatergic aberrations potentially secondary to dysfunctional gamma-aminobutyric-acid (GABA) receptors.6) Thus, several recent lines of evidence support neurotransmitter network abnormalities involving glutamate, GABA and dopamine.6) Interestingly, all these three neurotransmitter signaling pathways play a critical role in neuroplasticity.7)

Recent evidences suggest that aberrations in neuroplasticity8) and synaptic connectivity9,10) that are found in schizophrenia could be contributed by deficient brain-derived neurotrophic factor (BDNF).7) Blood serum3,11–15) and plasma levels16–18) of BDNF has been reported to be reduced in schizophrenia by several studies. Significant decreases in BDNF level have been noticed in hippocampus of people with schizophrenia.19,20) Further, these abnormalities in BDNF have shown associations with various clinical features15,21,22) of schizophrenia as well as antipsychotic treatment.23) Aberration in BDNF signaling pathways directly or circuitously influences neurotransmitters like glutamate and GABA, both of which are well understood to be altered in schizophrenia.24–26) These evidences suggest a central role for BDNF and its effector pathways in schizophrenia pathogenesis.

BDNF is most widely distributed neurotrophin in central nervous system and known to play a critical role in synaptic plasticity, memory persistence and hippocampal function,27,28) in addition to growth, proliferation and cell survival.29) These functions of BDNF are mediated by activation of the high-affinity receptor—tropomycin-related kinase B (TrkB).30,31) TrkB facilitates several functions like axon guidance, synapse formation and plasticity32–34) via activation of different pathways namely TrkB/phosphatidylinositol-3 kinase (PI3K), mitogen-activated protein kinase, phospholipase C gamma (PLC-γ), N-methyl-D-aspartate receptor (NMDAR) pathways.35) Relevantly, reduced and truncated TrkB protein has been observed in schizophrenia.36,37) The key components of these pathways are summarized below.

PLC-γ pathway

PLC-γ is significant for several cellular functions and implicated in brain disorders.38) TrkB promotes PLC-γ to release diacylglycerol (DAG) and inositol 1,4,5-trisphosphate (IP3) and in turn, IP3 is involved in differentiation of neuronal cells.39,40) Impairment of this pathway might be relevant for schizophrenia pathogenesis.41)

AKT pathway

Another pathway that is activated by BDNF corresponds to its neuro-protective function by activating AKT pathway.42) AKT inhibits glycogen synthase kinase 3 (GSK3)-β by phosphorylating it at SER-9 position and thereby promotes β-catenin towards division phase.43) This shows that GSK3-β is directly regulated by BDNF and could be affected in schizophrenia.

GABA synthesis and transportation

During regulation of GSK3-β, activated β-catenin induces GABA synthesis via activation of glutamic acid decarboxylase (GAD) 67.44,45) This indicates the indirect involvement of BDNF in regulation of GABA synthesis.46) GABA transporter 1 (GAT1) is involved in transportation of GABA outside the cell.47) GAT1 is modulated by BDNF through an intermediate pathway via PLC-γ.48) Impaired GAT1 has been reported in schizophrenia.38) Other factor conjoining GABA and BDNF is GAD67 which is also known to be altered in schizophrenia.49) Moreover TrkB is correlated to GAD67, as the decrease in TrkB is associated with decreases in GAD67 function.50)

Glutamate pathway

Synthesis of glutamate is chemical conversion of glutamine in the presence of glutaminase.51) GABA synthesis occurs by the conversion of glutamate to GABA in the presence of GAD67. VGLUT is glutamate transporter aids in glutamate release from glutamatergic cell and is influenced by BDNF. The pathway involved in the regulation of these processes is PLC-γ activation through TrkB.52)

Modeling BDNF can potentially facilitate comprehensive understanding of components of its signaling network and its contribution to molecular pathology of schizophrenia. Several previous attempts to understand the contribution of BDNF aberrations in schizophrenia have been conducted using animal models53–55) including drug treatment effects,56) role of microRNA and transcription factors57) etc. However, the underlying mechanisms of interactions are non-linear and the emergent biochemical function of a signaling network is often non-intuitive. Computational system biology provides an approach to study the frame of such a system and particularly help to understand the resulting behavior58). Further, biochemical modeling generates hypotheses which can be put to test experimentally in future studies. Moreover, in silico modeling approach provides additional advantages over direct experimentation on human or animal subjects59) such as being free from ethical considerations and providing an opportunity to perturb multiple parameters simultaneously to study their effect on the simulated environment. Hence, in this study we chose to apply in silico approach as opposed to in vivo/in vitro experimentation—although such an computational model can only provide restricted yet useful representation of reality.

We are unaware of any comprehensive quantitative modeling studies that have been published to elucidate the effects of BDNF on neurotransmitters like glutamate and GABA that are relevant for schizophrenia pathogenesis. Objective of current study is to construct a quantitative in silico model of signaling network integrating above mentioned pathways and consequently attempt to comprehend effect of BDNF and associated downstream signals on neurotransmitters. We applied biochemical modeling with a deterministic approach using knowledge from contemporary literature and parameters based on publicly available databases to understand links of BDNF with GABA and glutamate along with a set of important signaling factors that are known to be altered in schizophrenia (Fig. 1).

Fig. 1

Important components of the brain-derived neurotrophic factor (BDNF) network. Major conceptual components are shown. The nodes colored in red represent entities reported to be altered in schizophrenia. Edges represent interactions between the nodes and flow of signals, initiated by BDNF.

TrkB, tropomycin-related kinase; AKT, protein kinase B; NMDAR, N-methyl-D-aspartate receptor; PLC, phoshpholipase C; GABA, gamma-aminobutyric acid; CAM-II, Ca2+/calmodulin-dependent protein kinase II; CREB, cAMP response element-binding protein.

METHODS

Methods of current study were divided into three major parts: a) model construction for better understanding and representation of biological networks; b) simulation for estimating functionality of system in time domain; and c) analysis for obtaining results and converting it into comprehensible plots and matrices (Fig. 2). This pipeline consist of defining the model, initializing values, describing kinetic equations,60–71) verifying generated mathematical equations, outlining simulation parameters and algorithms as well as employing multiple analysis methods. These steps are explained sections described below.

Fig. 2

Workflow of the study methodology. (A) Multiple databases like Panther pathway, Reactome and BioModels were used for building model sub-components. (B) System Biology Markup Language (SBML) squeezer was used to generate the kinetic laws. Initial values were defined from BioGPS database. (C) The SBML model was imported in COPASI. Global constants were defined from BRENDA and DOQCS. Unavailable values were optimized after assigning parameters and defining ordinary differential equations. Analyses included steady state stability, perturbation time-course simulation and sensitivity analysis.

Model Construction

Representation of the network was built which aids in interpretation and analysing biological networks.72) Model was structured using System Biology Markup Language (SBML, a XML based language)73) in CellDesigner software (http://www.celldesigner.org/). SBML is a machine readable format for representing bio-models which can be simulated and analysed. CellDesigner supports all SBML versions and provide a user interactive interface to construct the model.74,75) CellDesigner also enables the import of SBML files from various sources and facilitates the export to various software.76) Publicly available databases such as Panther pathway, Reactome and BioModels were used for SBML model sub-components. Similarly, Uniprot, KEGG and BioSystems were used to understand the interactions as well as the functions of nodes involved.77) Various components of models were defined separately and merged to complete the interaction. For merging model or reactions already present in the database ‘merge model’ plugin was used. This plugin helps to merge SBML model manually as well as automatically.

Manual assignment of kinetic equation is cumbersome and highly error prone process, this was taken care by the plugin SBML squeezer.78) It is based on stoichiometry, the participating species and the regulatory relations. After SBML squeezer was invoked, kinetic equations were assigned to the reactions. These equations were cross checked and assigned names. For running simulations, we chose COPASI79) over CellDesigner as it provides better user interface and features to select a desired algorithm. Initial value, kinetic equations and parameters that were used in the study are provided as supplementary materials (Supplementary Tables 1–3). Concentrations of species in the network were initialized based on gene expression values for temporal lobe from BioGPS database.80)

Simulation

The SBML model that was created in CellDesigner was exported to the COPASI. This software was used for defining, simulating and analyzing computational models of biochemical network.79) Exported model was edited and updated in graphical user interface of COPASI.58) Ordinary differential equations (ODE) were generated in COPASI, details of which are provided as supplementary material (Supplementary Fig. 1). Values of global constants in network were set based on BRENDA and DOQCS.

Optimization of unknown parameters was done for quantities that were not available in databases (Supplementary Table 3). For optimization process combination of two algorithms was used—’random search’ for 10,000 iterations followed by ‘genetic algorithm’ for 300 iterations.81) We optimized each reaction sequentially. Depending upon the system configuration and number of parameters in the reaction for each optimization, run varied from 12 hours to 14 hours. The optimization task was set to time course and constraints boundaries were set between 0 to 100.

Time course simulation: Duration for simulation was 86,400 seconds (24 hours), with the interval period of 0.02 second. Comparison between the time points and the concentration of the species was performed. The method used was deterministic Livermore Solver for Ordinary Differential Equations with default parameter settings.

Analysis

Perturbation analysis was performed by varying the initial concentration of BDNF and recording the perturbed effect on the system. Four concentration levels were considered, which were 100%, 75%, 50% and 10% of initial concentration defined based on gene expression database BioGPS.

Sensitivity analysis was carried out by setting the task as time course with the ‘effect’ as non-constant concentrations of species and the ‘cause’ as initial concentration. Further, a ‘steady state’ analysis was performed which aims to deduce the stability of the system by analyzing behavior when the system achieves equilibrium. For this analysis parameters were set to default (newton=1, integration=1, back integration=1, accept negative concentration= 0) with 1,000 iterations. Jacobian matrix was calculated on complete system as well as reduced system. Steady state analysis also calculates the flux at equilibrium.

Graphical representation of results from the above analysis were generated using R programming software82) with integrated development environment—R Studio.83) The packages used were ‘plyr’ version 1.8.184) and ‘ggplot2’ version 1.0.1.85)

RESULTS

This study modeled key components of BDNF signaling network including as follows.

Calcium influx: Ca2+ influx and NMDAR activation by glutamate and BDNF; activation of cAMP response element-binding protein (CREB) pathway by triggering Calmodulin (CaM).

Protein kinase B/AKT pathway: Regulates apoptosis, stimulates cell cycle regulatory molecules and promotes GABA synthesis via AKT pathway.

Glutamate formation: Formation via glutminase and activation of glutamate transport via BDNF-TrkB mediated PLC-γ pathway.

GABA formation: Conversion of glutamate to GABA via β-catenin and GAD67; transportation of GABA via BDNF-TrkB interaction with VGAT (Fig. 3).

Fig. 3

Network fabricated in CellDesigner. Components of the modeled network were (A) calcium influx: Ca2+ influx and NMDAR activation by glutamate and brain-derived neurotrophic factor (BDNF); activation of CREB pathway by triggering calmodulin. (B) Protein kinase B/AKT pathway: Regulates apoptosis, stimulates cell cycle regulatory molecules and promotes gamma-aminobutyric acid (GABA) synthesis via AKT pathway. (C) Glutamate formation: Formation via glutminase and activation of glutamate transport via BDNF-tropomycin-related kinase B (TrkB) mediated PLC-γ pathway. (D) GABA formation: Conversion of glutamate to GABA via β-catenin and glutamic acid decarboxylase (GAD) 67; transportation of GABA via BDNF-TrkB interaction with VGAT.

NMDAR, N-methyl-D-aspartate receptor; CREB, cAMP response element-binding protein; AKT, protein kinase B; PLC, phospholipase C; VGAT, vesicular GABA transporter.

Short term effect of BDNF was observed on GABA (Fig. 4). Progression in activation of species were analysed. Peak of activation of BDNF was achieved at ~80 seconds. Time course of the system stabilized after ~100 seconds. BDNF and TrkB interaction produced a catalytic effect to various reactions including conversion of glutamate to GABA via GAD67, transportation of glutamate and GABA, activation of NMDAR receptor and activation of AKT. The intracellular GABA concentration peaked at ~0.1 seconds and peak of extracellular GABA concentration was observed at ~10 seconds. Peak in concentration of activated AKT was observed in ~3 seconds and plateaued at ~7 seconds. Maximum of GABA concentration was obtained at ~50 seconds inside GABAergic cell. Most of the GABA produced was secreted outside the cell and only 5% of the total GABA remained inside the cell. Beta-catenin activation was rapid which initiated conversion of GABA.

Fig. 4

Time course plot for gamma-aminobutyric acid (GABA). Short term effect of brain-derived neurotrophic factor (BDNF) on GABA: time course of 100 seconds with 5,000 intervals of 0.02 seconds depicted on log scale. *Activated; red arrow, inhibition; green arrow, activation.

TrkB, tropomycin-related kinase B; AKT, protein kinase B; GSK3-β, glycogen synthase kinase 3; CTNNB1, catenin beta-1; GAD67, glutamic acid decarboxylase 67.

Short term effect of BDNF on glutamate was analyzed and plotted (Fig. 5). Time course depicts formation as well as release of glutamate from the cell. Activated NMDAR was ~80% of the total volume of NMDAR. Active PLC-γ plateaued at ~5 seconds. VGLUT was activated within ~1 seconds of BDNF activation that simultaneously increased extracellular glutamate leading to depletion in concentration of intracellular glutamate.

Fig. 5

Time course plot for gluamate. Short term effect of brain-derived neurotrophic factor (BDNF) on glutamate: time course of 100 seconds with 5,000 intervals of 0.02 seconds depicted on log scale. *Activated; green arrow, activation.

TrkB, tropomycin-related kinase B; PLC, phopholipase C; VGLUT, vesicular glutamate transporter.

Perturbation analysis was done by varying BDNF concentration gradually from 100% to 10% (Fig. 6). Perturbed level of BDNF, i.e., 100%, 75%, 50% and 10% of initial concentration influenced downstream signaling processes. Lower BDNF concentration showed delay in several processes. The time at which maximum activated concentration of PLC-γ was observed in the presence of 100% BDNF (~3 seconds) was reduced ~4 fold with lowest level of BDNF. Presence of 10% BDNF level disrupted the activation of AKT and required ~11 seconds more for stabilization. Rate of glutamate conversion to GABA and glutamate transport was decreased at lowest BDNF level. Time required for transport of glutamate from glutamatergic cell to extracellular space was almost doubled when BDNF was reduced to 10%. Rate of formation of GABA was slower with lower levels of BDNF. Concentration of activated NMDAR obtained within the time at which maximum concentration was observed in the presence of 100% BDNF (0.6 seconds) reduced ~2.5 fold with lowest level of BDNF.

Fig. 6

Effects of brain-derived neurotrophic factor (BDNF) perturbations. Time course plots on log-linear scale at varying BDNF concentrations for (A) PLC-γ, (B) AKT, (C) glutamate in GABAergic cell, (D) glutamate in glutamatergic cell, (E) GABA extracellular, and (F) AKT for 100 seconds with 5,000 intervals of 0.02 seconds. Colors of the lines from dark to light pink correspond to the percentage of initial BDNF concentration with respect to database derived expression levels in temporal lobe, i.e., 100%, 75%, 50%, and 10% respectively.

PLC, phospholipase C; AKT, protein kinase B; GABA, gamma-aminobutyric acid; NMDAR, N-methyl-D-aspartate receptor.

Sensitivity analysis was performed on the model to identify influential species (Supplementary Table 4). Matrix obtained from analysis shows association between the species and depicts the strength of the interaction between them (Supplementary Fig. 2). Among the stronger interactions that were found, the relationship between CREB and CREB-regulated transcription coactivator 1 (CRTC1) was strongest followed by TrkB-BDNF and calcium-CaM. Calcium concentration was more inversely related to concentration of PLC-γ than that of CaM. Sensitivity calculation showed that in addition to TrkB and BDNF, CaM, CREB and CRTC1 have significant impact on the system (Fig. 7).

Fig. 7

Sensitivity analysis. Bar plot of the sensitivity analysis: x axis displays the species and y axis shows sensitivity. In addition to BDNF-TrkB, calcium, calmodulin, CRTC1 and CREB show high sensitivity.

AKT, protein kinase B; BDNF, brain-derived neurotrophic factor; CAM2, Ca2+/calmodulin-dependent protein kinase II; CREB, cAMP response element-binding protein; CRTC1, CREB-regulated transcription coactivator 1; CTNNB1, catenin beta-1; GAD67, glutamic acid decarboxylase 67; GSK3, glycogen synthase kinase 3; NMDAR, N-methyl-D-aspartate receptor; PIP2, phosphatidylinositol 4,5-bisphosphate; PLC, phospholipase C; PPP3, protein phosphatase 3; TrkB, tropomycin related kinase B; vGAT, vesicular GABA transporter; VGLUT, vesiccular glutamate transporter.

We performed kinetic stability analysis—a linear steady state stability analysis based on the eigenvalues. An equilibrium steady state with near-zero fluxes could be attained. All negative values in reduced Jacobian showed that steady state was asymptotically stable. Twenty-four real eigen values were detected with largest real part being −0.07. This indicates that these species can preserve same orientation/vector in the model.

DISCUSSION

The study provides an in silico quantitative model of biochemical network constituting the key signaling molecules implicated in schizophrenia pathogenesis. Focus of the study was to analyze aspects of BDNF function on neurotransmitters specifically GABA and glutamate based on publicly available data or existing literature. We present an integrative scaffold for the processes with a deterministic approach.

BDNF could be one of the candidate genes for schizophrenia8) because of its direct involvement in neurodevelopment,46) neuroplasticity,7) interaction with neurotransmitters86,87) and its association with clinical features such as age of onset,22) clinical subtypes,15) symptom severity21,88) and drug response23)—both at genetic89,90) as well as blood level.11,17) However, these clinical associations have not always been consistent.91,92) Interestingly, antipsychotics can affect BDNF levels93,94) potentially mediated through epigenetic mechanisms.89,95) More studies are required to clearly elucidate the role of BDNF and its effector pathways in pathogenesis of schizophrenia.

Our model demonstrated that shortly after activation of BDNF, it triggers other enzymes to carry out processes such as transcription through CREB phosphorylation which promotes cell division, promotion of anti-apoptotic pathway through AKT42) which in turn triggers glutamate conversion to GABA44) (by inhibiting GSK-β and activating β-catenin), transport of GABA via vGAT, as well as direct and indirect transport of glutamate via VGLUT and activation of PLC-γ48) respectively. Reduction of BDNF implicated delay in rate of PLC-γ, NMDAR and AKT activation, GABA formation, glutamate and GABA transport. These key players could potentially influence brain signaling and neuroplastic function via temporal modulation of time sensitive neurotransmitter release. However, as observed in our simulations, stable concentrations could be achieved back as system tries to neutralize the perturbation effect through compensatory interactions.96) Hence, the functional impact of these perturbations would be dependent upon the time domain of its availability for downstream processes. Other players with high impact on system were calcium, CaM, CREB and CRTC1.

It should be noted that BDNF is only one of several signaling molecules that participate in pathophysiology of this disorder. Isolated attention to BDNF and its effects could be considered as a potential limitation of the study and future studies should expand the scope to dopamine,4) serotonin97) and other neurotransmitters.98) Further, analysis of the neural system in schizophrenia is challenging due to involvement of vast, varied and substantially complex mechanisms.7,25,99) Critical regulatory mechanisms might be omitted within the scope of a single study due to methodological limitations which might lead to inaccurate or incomplete model representations. For instance, the influence of BDNF on dopamine via D2-like/D3-like receptors100) has not been modeled in this study. Moreover, there are schizophrenia risk genes other than ones modeled in this study (i.e., AKT, BDNF, and NMDAR) such as DISC-1, neuregulin/ErbB4, dysbindin etc. which are involved in regulating neuroplasticity.7) With fast-growing computational advancements, future studies can possibly model genetic and transcriptional activity to characterize specific bimolecular mechanisms that drive neural changes in this disorder. However, building quantitative models of such extensive scope of species with multiple compartments poses significant methodological challenges due to missing links in various aspects of current knowledge i.e. clear elucidation of reactions in pathways and their molecular interactions, transport and compartmentalization of species, specific kinetic parameters, concentration/expression values for various cell types etc.

In relevance to schizophrenia pathogenesis, the key insights from current study are increased glutamate in hippocampus observed in this disorder24) could be contributed by an indirect pathway wherein reduced BDNF suppresses GAD67-mediated GABA synthesis. Further, deficient BDNF corresponds to impaired transport via VGLUT, thereby further increasing the intracellular glutamate in glutamatergic cells. BDNF also alters the calcium dependent neuroplasticity via NMDAR modulation. In summary, this study provides mechanistic insights into putative contribution of deficient BNDF towards alterations in neurotransmitters and neuroplasticity that are consistent with current understanding of the schizophrenia. Future works need to evaluate the relationship of BDNF perturbations with dopaminergic aberrations and also, extend the scope to include other neuroplasticity pathways regulated by schizophrenia risk genes such as DISC-1, neuregulin/ErbB4, dysbindin etc.