Modeling the activity of the dopamine signaling pathway by combination of analog electrical circuit and mathematical approaches.

This paper demonstrates the application of the system biology principles on the example of the dopamine signaling pathway in neurons. Presented model is based on two approaches - cytomorphic electronic circuits and mathematical modeling. Transcription and phosphorylation of DARPP-32 was modeled by analog circuit, based on well-known approaches presented in [1]. It was shown that application of circuit helps to receive signal oscillations that close to described ones in real biological systems. This combination on the one hand gives possibility to simplify calculations, on another to show this signaling pathway dynamics. The expected effect of changes in the functioning of calcium channels is considered, and the mathematical model of the interaction of system components is proposed. The average frequency of calcium current oscillations due to the presence of dopamine was 30 Hz in presented model, that is consistent with the literature, where the frequency of such oscillations is up to several tens of Hz. All presented results shows good correlation with known data, which already published today.

Introduction

To study complex biological systems, scientists often turn to special subdivisions of modern biology — system and computational biology and biochemistry, which use the advances of computer science, physics, and apply mathematics to solve problems [2].

Today, there are many approaches to the description of cellular processes at the genetic level, which allows us to extend the idea of considering a single biological system as a complex of several subsystems [3]. For example, the analogy between synthetic biological and electrical circuits is actively discussed today [1]. The processes of electrical circuit's elements operation simulate natural processes of regulation in the cell. An analogy was also shown between the dynamics of biochemical reactions and the flow of electrons in transistor, which makes it possible to expand the field of application of the idea of using electronic circuits with the help of transistors. Analog electronic flow in subthreshold transistors and analog molecular flux in chemical reactions obey Boltzmann exponential laws of thermodynamics and are described by astoundingly similar logarithmic electrochemical potentials. Therefore, cytomorphic circuits can help to map circuit designs between electronic and biochemical domains. Thus, in connection with many analogies between different processes in cells and transistors [1], a circuit implementation of dynamic processes of transcription and translation in a living cell was proposed.

Signaling networks play the central role in the regulation of processes in a single cell and in the entire body. The dopaminergic signaling pathway was chosen as the object of the study (Figure 1). Dopamine is responsible for the occurrence of certain cellular responses, for example, a change in the conductivity of calcium channels in neurons [4]. It also plays an important role in ensuring cognitive activity, and the lack of dopaminergic transmission is one of the causes of various cognitive diseases, such as Parkinson disease. Symptoms include motor tremor, difficulty initiating actions and eventually catatonia and anhedonia. On the other hand, when the dopamine neurons behave normally, they provide brief bursts of dopamine to the neocortex and other brain areas when a reward (unexpected pleasure) occurs and a diminution of activity when less than expected reward is experienced.

Figure 1

Dopaminergic signaling pathway. DARPP-32 phosphorylation can be controlled by synaptic release of dopamine, and phosphoprotein is involved in the regulation and management of dopamine-induced trans-synaptic effects, which is achieved by the association of D1 receptors with the activation of adenylyl cyclase [8], triggering cAMP synthesis [9, 10]. Under normal conditions, DARPP-32 is phosphorylated by the residues Thr34 by PKA (protein kinase A [12]), Thr75 by CDK-5 (Cyclin-dependent kinase 5) [13,14], Ser97 and Ser130 [15]. Dopamine activation of the D1 receptor and, accordingly, of this signaling pathway increases the level of P-DARPP-32-Thr34, which is a PP-1 (protein phosphatase 1 [27]) inhibitor and stimulating the PKA-sensitive regulatory subunit protein phosphatase 2A, which, in its turn, dephosphorylates Thr75 residues. Dephosphorization of the Thr75 residue directly leads to inhibition of PKA. Thus, in the dopaminergic signaling pathway there is a cycle with positive feedback. One of the effects of dopamine exposure and, accordingly, activation of the signaling pathway through D1 receptors is a change in the properties of ion channels, for example, their conductivity [17].

Functioning of DARPP-32

DARPP-32 phosphoprotein (dopamine and cAMP-regulated phosphoprotein, Mr 32,000 kD) plays an important biological role in dopamine-sensitive neurons [5, 6]. Recent combination of experimental and modeling studies suggests that DARPP-32 is a robust integrator of signaling [30] whose main role may be to increase the reliability in decoding the information mediated by dopamine, as well as other inputs (such as glutamate) [6]. It is localized mostly to regions that receive dopaminergic innervation is and present in virtually all medium-size spiny neurons of the striatum. DARPP-32 plays a major role in the dopaminergic signaling pathway. Regulation of the state of DARPP-32 phosphorylation provides a mechanism for integrating information arriving at dopaminoceptive neurons, in multiple brain regions, via a variety of neurotransmitters, neuromodulators, neuropeptides, and steroid hormones [7]. DARPP-32 phosphorylation can be controlled by synaptic release of dopamine, and phosphoprotein is involved in the regulation and management of dopamine-induced trans-synaptic effects, which is achieved by the association of D1 receptors [7] with the activation of adenylyl cyclase [8], triggering cAMP synthesis [9, 10]. cAMP in turn activates the protein kinase A [10] through release of the catalytic subunits when levels of cAMP rise in response of dopamine signal. Catalytic activity of adenylyl cyclase, and thus effectiveness of PKA phosphorylation [11], is regulated by different molecules, including G-proteins, capable of both stimulating and inhibiting its work. Under normal conditions, activity of DARPP-32 depends on the state of phosphorylation at multiple regulatory sites, including Thr34 [12], Thr75 [13,14] and additionally Ser97, Ser130 [15]. Phosphorylation pattern itself depends on the dynamic balance between activation of phosphatases and protein kinases. Phosphorylation at Thr34 by PKA converts DARPP-32 into an inhibitor of Protein Phosphatase-1 (PP-1), phosphorylation at Thr75 by cyclin-dependent kinase 5 converts DARPP-32 in an inhibitor of PKA. It turns out that the mechanism of influence of CDK-5 on the dopamine signaling pathway in neurons is strongly related to the work of the coronin 1 protein, which interacts and binds to the G protein, which in turn directly affects the synthesis of cAMP via adenylyl cyclase. Dopamine activation of the D1 receptor and, accordingly, of this signaling pathway increases the level of P-DARPP-32-Thr34, which is a PP-1 [27] inhibitor, and stimulating the PKA-sensitive regulatory subunit protein phosphatase 2A [16], which, in its turn, dephosphorylates Thr75 residues. Dephosphorization of the Thr75 residue directly leads to inhibition of PKA. Thus, in the dopaminergic signaling pathway there is a cycle with positive feedback (Figure 1). One of the effects of dopamine exposure and, accordingly, activation of the signaling pathway through D1 receptors is a change in the properties of ion channels, for example, their conductivity [17].

As for existing models of DARPP-32 activity, approach is mostly the same - all chemical reactions are described as protein-protein interactions or as enzymatic reactions written in Michaelis-Menten form [22]. Solutions of corresponding first-order differential equations can be compared to the experimental data [9]. In [22] robust simulation of DARPP-32 signaling pathway was obtained with use of reaction rate constants from the scientific literature. In [9] more detailed dopamine impact analysis was provided. Simulations were used to compare the effect of different dopamine receptors in neurons.

In this paper, we present the application of the hybrid approach for modeling of intracellular signaling pathways based on the combination of analog electrical circuit and mathematical models. Such an approach allows one to obtain a reliable mechanism for modeling metabolic pathways based on multiple analogies of the thermodynamics of processes at MOSFET and in biochemical reactions [1]. In this work we will be focused on the dynamics of calcium channel current behavior when the concentration of dopamine is changing during the signaling cascade induction by the activation of dopamine receptors of type D1 is considered.

Materials and methods

Analog circuit for modeling transcription and translation of DARPP-32

The use of transistors circuit allows us to determine the transfer function of the entire circuit. It also helps to consider the influence of molecular noise [18], bringing our model closer to the real dynamics of changes in concentration in a living cell [19]. Also at the moderate precision of computation seen in cells, analog computation is significantly more efficient in its use of energy, time and space (molecular count and part count) than digital computation. Therefore, analog and collective analog computation is likely to be more practical and scalable for synthetic biology than purely digital computation [1]. The input parameters are the activator and repressor voltages, the difference of which considers the activity of the transcription factor of the particular protein. The voltage, which as a result of transformations is obtained at the output, is proportional to the concentration of the translated protein. This task also implies the solution of the Michaelis-Menten equations, i.e. ODE of the 1st and 2nd order. Electrical circuits that solve such equations are already developed and are described in detail [19]. The interest of this approach consists in obtaining a reliable method for predicting a biological response to changes in one or another component of a complex system without interfering with the cell itself. It is also assumed to take into account molecular noise, which is difficult to achieve in other methods [1, 20, 24].

According to the analogy of the processes of transcription and translation [1], the circuit reproducing the dynamics of translation of the phosphoprotein DARPP-32 was designed (Figure 2). Iind, Isup are the currents proportional to the concentration of activators and repressors of the transcription of the protein of interest, Xt is the total concentration of the transcription factor, Dt is the maximum concentration of the protein, Kd is the dissociation constant.

Figure 2

Scheme that reproduces the dynamics of the DARPP-32 [1]. Iind, Isup are the currents proportional to the concentration of activators and repressors of the transcription of the protein of interest, Xt is the total concentration of the transcription factor, Dt is the maximum concentration of the protein, Kd is the dissociation constant.

Mathematical model

One of the effects of dopamine exposure and, accordingly, activation of the signaling pathway through D1 receptors is a change in the properties of ion channels, for example, their conductivity. So, PKA, which work is modulated by DARPP-32, enhances the work of the Cav1.2 channels at high concentrations of Ca2+, weakening the connection of the calmodulin with the terminal tail of the C1 channel [17]. DARPP-32 and its effect on the dopaminergic signaling pathway (phosphorylation of DARPP-32 using CDK-5 and further inhibition of PKA) are responsible for the change in the conduction of calcium channels in neurons, and therefore affect the signal transmission in the brain. This conclusion was made by scientists back in 1999 [13], by analyzing these channels in individual striatal neurons using the patch-clamp technique.

In this paper, we are presenting a mathematical model of the calcium channel when changing the concentration of dopamine. The chemical reactions between the components of the signal path described above were calculated in accordance with the Michaelis-Menten kinetic equations by analyzing and numerically solving differential Eq. (1). Considering dopaminergic signaling pathway we have several systems of Mehaelis-Menten equations (E – PKA, CDK-5; S - DARPP-32; P – P-DARPP-32-Thr34, P-DARPP-32-Thr75 respectively). The necessary rate constants for the direct and reverse reactions of this path were taken from the corresponding articles [21].E is an enzyme, S is a substrate, ES is an enzyme-substrate complex, P is a reaction product, k1 is a reaction rate constant for the formation of an enzyme-substrate complex from an enzyme and substrate, k-1 is a reaction rate constant for the dissociation of an enzyme-substrate complex into an enzyme and substrate, k2 — rate constant for the transformation of an enzyme-substrate complex into an enzyme and product.

The main parameters determining the behavior of the control system of interest are the concentrations of PKA, CDK-5, DARPP-32 and its phosphorylated forms on threonine residues 34 and 75, dopamine: phosphoprotein DARPP-32 as substrate S, kinases PKA and CDK-5 as enzymes E, P-DARPP-32-Thr34 and P-DARPP-32-Thr75 as products of reaction P in Eq. (1). As a first approximation, the concentration of PKA and CDK-5 in the cell is considered to be much higher than the concentration of the other components of the system, which simplifies the system of Mehaelis-Menten differential equation and allows you to write an approximate analytical solution. The deviation of the model with such an assumption from the models with numerically solved DE [9, 22] is insignificant. In addition, in this model, using the methods of system biology, transcription and translation of components of the system are taken into account.

In (2) A and R stand for activator and repressor concentration respectively, δ is coefficient of concentration decay, f1,2 are Hill functions in which K1 and K2 are the maximal expression rates when the activator protein is in excess, KA and KR are the expressions in the absence of activator or repressor, γ1, γ2 and γ3 are coefficients related to the affinity between the proteins and the promoter regions, and n is the Hill coefficient [23].

To calculate the concentration of DARPP-32 produced in a cell at some point over time, it is necessary to understand how the concentration of activators and repressors of a particular protein transcription changes in the cell. According to [23], their concentrations are determined by solutions of the differential Eq. (2). Also one can analyze these equations more thoroughly, paying particular attention to the existence of a solution, the features of its trajectories and its stability. With certain initial conditions [23] we get stable trajectories (Figure 3). Detailed explanation and discussion is presented in [24]. Considering dopaminergic signaling pathway, A and R are activator and repressor concentration of DARPP-32, which in turn acts as substrate S in Eq. (1). So dynamics of Eq. (2) indirectly affects products P of Eq. (1).

Figure 3

Changes in the concentration of DARPP-32 activator and repressor in the cell over time with respective phase trajectories of system (2) for different repressor decay parameter δr (A: δr = 0.5, B: δr = 1.0, C: δr = 1.08, D: δr = 1.5).

According to [24], resulting circuit voltage, which is proportional to the concentration of the translated DARPP-32, depends on (3):

C is a constant depending on the external conditions of the system, is the difference of the voltages applied to the circuit (Figure 2) corresponding to the voltages of the activator and repressor. Circuit input currents (Figure 2) are generated respectively based on solutions of Eq. (2) (A and R), thus in circuit experiment we obtain concentration of DARPP-32 based on its activator and repressor dynamics and experimental modeling of the protein translation process using electrical circuits.

In this model, the effect of changing the PKA and PP-1 concentration as a result of interaction with DARPP-32 is to change the conductivity of the calcium channel in accordance with the Hodgkin-Huxley model, according to which the calcium current is defined as:

- membrane permeability, V - membrane potential, E - membrane resting potential, m = m(V) and h = h(V) - coefficients responsible for channel activation and inactivation, non-linearly dependent on V. Expressions for them were taken from [25]. Knowledge of changes in permeability of the channel due to changes in concentration of PKA and PP1 gives us the relationship between changes in the concentration of DARPP-32 and changes in calcium current.

So changing in DARPP-32 concentration will affect permeability of the calcium channel introduced in Eq. (4). This equation was solved in conjunction with the control for the concentration of calcium [31]:

[Ca]i is the concentration of intracellular calcium [Ca]∞ is the equilibrium concentration of intracellular calcium.

Parameters and initial conditions of variables of mathematical model

All rate constants for biochemical equations and initial conditions for substrate, enzyme and product concentration were taken from data, based both on mathematical simulations [26] and biological experiments [9, 22]. Those initial conditions of variables can be found in supplementary materials of our paper in corresponding section. For convenience, initial conditions of variables used in this work can be found in Table 1.

Table 1

Initial conditions for variables of mathematical model. Additional data can be found in supplementary materials of this paper.

Variable	Value	Unit	Reference	Comment
Da	10	μM	[22]	Dopamine
DARPP-32	50000	nM	[9]	Dopamine and cAMP-regulated phosphoprotein, Mr 32,000 kD
P-DARPP-32-Thr-34	400	nM	[9]	DARPP-32, phosphorylated by the residue Thr34 by PKA
P-DARPP-32-Thr-75	12000	nM	[9]	DARPP-32, phosphorylated by the residue Thr75 by PKA
PKA	1,2	μM	[22]	Protein Kinase A
CDK-5	1,8	μM	[22]	Cyclin-dependent kinase 5
PP1	5	μM	[22]	Protein Phosphatase-1
cAMP	5000	μM	[22]	Cyclic adenosine monophosphate

Results

System (2) solutions for different δr can be either periodic (0.5, 1) or damped (1.5), as shown in Figure 3. The frequency of oscillations of the activator and repressor concentration depends on the repressor decay parameter δr, decreasing in time. Also noticeable is the change in the nature of the concentration dependence of this parameter - with its increase, the periodic dependence is replaced by a decaying periodic, and then a decaying non-periodic one. The frequency for the case with δr = 0.5 is 0.09 Hz, for δr = 1.0 it is 0.14 Hz and for δr = 1.5 there are no oscillations.

Such behavior suggests the idea of studying the joint dynamics of the concentration of activator and repressor DARPP-32 over time. Figure 3 shows the trajectories of solutions in phase space, which indicate the stability solutions of system (2). The obtained data allow us to say that with an increase in the parameter δr, corresponding to the degradation of the repressor, the system tends to asymptotic stability. The critical δr value is approximately 1.08, solutions of system (2) and corresponding phase portrait is shown in Figure 3 as well. At smaller values of the parameter, the solutions of system (2) will be periodic, which corresponds to non-asymptotic stability. The profile of phase portrait evolution over time is presented on Figure 4, where critical value of repressor decay parameter corresponds to drastic change in external part of the profile's envelope.

Figure 4

Changes in the phase portrait of system (2) as the repressor decay parameter δr changes from 0.5 to 1.5.

Using Eqs. (3), (4), and (5) one can obtain dependencies of concentration of phosphorylated DARPP-32 over time. The resulting change in the concentration of DARPP-32 phosphorylated at residues of threonine 34 and 75 (Figure 5) corresponds to a numerical calculation based on data from real organisms reaching concentration values within 12 μM compared with 1–10 μM range of corresponding concentrations from literature [22]. The behavior of these dependencies (Figures 5B - 5D) reflects the same relation to repressor decay parameter – with increasing δr oscillations of phosphorylated protein concentration changes in similar way to change in activator and repressor concentration over time. In this work we compared dynamics of protein concentration and corresponding calcium current dynamics in case of the absence of taking into account the dependence of the concentration of the translated protein. In this case, the voltage difference in Eq. (2) is constant (Figure 5A). Еhe comparison shows that in the case of a more complete model, the amplitude of the oscillations of the calcium concentration in time are within the same limits as in the base model. Moreover, when the repressor decay parameter increases, the concentrations of phosphorylated proteins cease to oscillate and the dynamics of their behavior tends to the baseline dynamics.

Figure 5

Changes in the concentration of phosphorylated DARPP-32 over time (A: and for different repressor decay parameter δr (B: δr = 0.5, C: δr = 1.0, D: δr = 1.5)).

The dependences of the calcium current in the presence and absence of dopamine on time were obtained, as well as the joint dynamics of changes in the concentration of calcium, membrane potential and calcium current (Figure 6). The effect of intracellular calcium on the functioning of the PKA itself was neglected. As it was expected, the dynamics of calcium current tends to baseline one (Figure 6A) with increasing of repressor decay parameter (Figure 6D). The frequency of oscillation is similar, but amplitude is different. For intermediate in terms of δr case several distinct peaks in calcium current are observed.

Figure 6

The dependence of calcium channel current in the presence and absence of dopamine over time. () and for different repressor decay parameter (B: δr = 0.5, C: δr = 1.0, D: δr = 1.5)).

As a result of supplying the electric circuit (Figure 2) of the voltages Vrep and Vact, obtained by solving system (2) for different values of the system parameter δr (0.5, 1.0, 1.5). The voltage read from the circuit will be proportional to the concentration of phosphoprotein DARPP-32. In accordance with this, the dynamics of the change in the concentration of P-DARPP-32-Thr-34 and P-DARPP-32-Thr75 was calculated. Eqs. (4) and (5) were solved with the obtained data.

A comparison was made of the obtained data with the results of mathematical modeling considering the DARPP-32 translation using (3). Dependencies are built together (Figure 7). Average frequency of oscillations of calcium current due to the presence of dopamine is 30 Hz. Both frequency and amplitude is close in both cases, as well as dynamics of calcium oscillations in general, such as profile envelope and peak times (Figures 7A, 7B).

Figure 7

The dependence of calcium channel current over time obtained as a result of a complete mathematical modeling, where concentration of DARPP-32 is calculated using only Eqs. (1) and (2) (blue line) and combinations of modeling and experiment, where concentration of DARPP-32 is calculated using data from circuit and Eqs. (1) and (2) (red line) for different repressor decay parameter (A: δr = 0.5, B: δr = 1.0, C: δr = 1.5).

Discussion

To analyze the correlation of the dependencies obtained, the coefficient of determination (Table 2) was calculated. The obtained correlation coefficients show a significant relation between the calcium current dependencies obtained as a result of the experiment and from the simulation.

Table 2

The values of the coefficient of determination (R-square) for dependencies with different

δ_r	R2
0.5	0.72
1.0	0.73
1.5	0.96

Based on the calculated values of the coefficient of determination (Table 2), we can speak about the presence of a substantial relationship between the results obtained as a result of modeling and experimenting with an electric circuit.

A circuit simulating biological objects similar to those considered in this paper is quite simple to assemble, and their use in such problems is justified by multiple analogies of the thermodynamics of processes at MOSFET and in chemical reactions [1]. Moreover, a similar scheme can be extended further, including new components that allow describing more complex dynamics. Such studies will help to understand better the functioning of metabolic signaling pathways and corresponding associated diseases.

The profile and amplitude of the calcium current, constructed using formulas (4formulas (4) and (5)(5) is comparable to the data from literary sources [28]. Also, data on changes in the phosphoprotein DARPP-32 component in the presence or absence of dopamine are consistent with information presented at paper [22]. The average frequency of calcium current oscillations due to the presence of dopamine was 30 Hz, which is also consistent with the literature, where the frequency of such oscillations is up to several tens of Hz [29]. Such changing in profile with repressor decay parameter changing is expected - the bigger it gets; the less experiment results diverse from situation with constant concentrations of activator and repressor. It occurs because of neglecting repressor concentration change due to Eq. (2). The situation with oscillation of activator and repressor concentrations leading to corresponding changes in phosphorylated DARPP-32 concentration, shows partial decrease in calcium current oscillation intensity in neurons in response of dopamine signal, which is consistent with data, presented on the scientific literature [26]. Amplitude of calcium current oscillations reaches several pA, which is consistent with described in the scientific papers direct measurements data as well [28]. Thus, model gives possibility to investigate influence of activator and repressor concentrations dynamics to the DARPP-32 concentration and calcium channel current in the end with correlation to the dopamine concentration.

It should also be noted that the biosynthesis of proteins involved in signaling pathways has a significant effect on the frequency of Ca-channel oscillations. From the presented data, it can be seen that, according to the model, a change in the ratio of activator/repressor concentrations causes changes in the frequency of corresponding current oscillations. Further research will be aimed at studying the dynamics of biosynthesis of all proteins involved in this signaling pathway and assessing its effect on Ca channels' work.

Conclusion

This results overall show ability to relatively accurate simulate activity of such biological objects as signaling pathways using analog circuit approach via cytomorphic transistor-based circuits equivalents. This tool can allow scientists not only to model some biochemical reactions but also consider more realistic simulation including biological molecular noise and analyze the robust response of a complicated network of reactions to an input signal in this way.

Based on the good correlation with known experimental data, we can conclude that presented in paper model that based on combination of cytomorphic analog electronic circuit and mathematical approaches can show satisfactory results with good correlation with previously published in literature data.

Declarations

Author contribution statement

Shumilov A.V.: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Wrote the paper.

Gotovtsev P.M.: Conceived and designed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

This work was supported by the , Thematic Plan 1.13.

Data availability statement

Data included in article/supplementary material/referenced in article.

Declaration of interests statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.