Modeling domain formation of MARCKS and protein kinase C at cellular membranes

*Correspondence: sergio.alonso@ptb.de Physikalisch-Technische Bundesanstalt, Abbestrasse 2-12, 10587 Berlin, Germany Abstract Background: Phosphorylation and dephosphorylation of proteins are mechanisms of activation and deactivation which regulate many cellular processes. Both mechanisms have been usually described in well mixed environments. MARCKS is a protein which binds to the membrane by electrostatic interaction. It is translocated from the membrane and phosphorylated by Protein Kinase C. Back in the cytoplasm the translocated MARCKS proteins are dephosphorylated by the enzyme phosphatase and can reattach to the membrane. These three processes (membrane binding, translocation, dephosphorylation) give rise to a cyclic dynamics known as the myristoyl-electrostatic switch.


Background
Living cells process external perturbations by biochemical pathways which are regulated by complex metabolic networks. If a receptor at the membrane is activated by an external signal, a cascade of biochemical reactions and interactions inside the cell is triggered. Proteins, enzymes and small molecules interact by mutual activation and inhibition. Multiple molecules participate in metabolic networks, where phosphatases and kinases are particularly active in cascades that involve phosphorylation and dephosphorylation of proteins involved in the signaling [1].
The number of components of such networks is usually very large. To model the relations among the components, the interactions are simplified and the pathways are http://www.epjnonlinearbiomedphys.com/content/2/1/1 described by ordinary differential equations of the component concentrations in the cell. However, some molecules and reactions are restricted to particular compartments of the cell [2]. As a result, the spatio-temporal distribution of proteins in the cell is important [3]. The communication between steps of the pathways restricted to different compartments is typically achieved by diffusion and active transport of the involved molecules. For sufficiently high molecular concentrations, the spatial distribution can be modelled by reaction-diffusion equations. These equations permit also the implementation of spatially inhomogeneous conditions, which reflect the differences among separate compartments inside the cell.
The myristoylated alanine-rich C kinase substrate (MARCKS) is a charged and unfolded protein which appears in living cells at high concentrations of around 10μM [4]. This protein is active in many processes occurring at the cell membrane demanding the regulation of cytosketetal dynamics, e.g. cell movement [5], phagocytosis [6], exocytosis [7], and chemotaxis [8]. MARCKS is a hydrophilic protein which can bind to the membrane by electrostatic interactions with acidic lipids like phosphatidylinositol 4,5-biphosphate (PIP 2 ) or phosphatidylserine (PS) [4]. The functions of MARCKS are correlated with its ability to sequester PIP 2 at the membranes and to avoid thereby its hydrolysis by phospholipase C (PLC). Thus, MARCKS prohibits the generation of two important second messengers: inositol 1,4,5-triphosphate (IP 3 ) and diacylglycerol (DAG), which participate in the regulation of calcium (Ca 2+ ) in the cell.
MARCKS proteins lose their affinity to membranes when they are phosphorylated by protein kinase C (PKC). The phosphates reduce the positive charge of the protein and cause the unbinding from the membrane to the cytoplasm. In the cytoplasm, phosphatases remove the phosphates again from the protein. Consequently, MARCKS can bind again at the membrane. This cyclic dynamics of MARCKS binding and unbinding is known as the myristoyl-electrostatic switch [9]. It is controlled by Ca 2+ which activates and promotes the binding of PKC at the membrane, where it phosphorylates MAR-CKS proteins. Therefore, phosphorylation and dephosphorylation of MARCKS appear in different regions of the cell, in contrast to the earlier Goldbeter-Koshland model of protein phosphorylation and dephosphorylation [10,11], where these processes are assumed to occur in a well-mixed environment inside the cell. The localization of both processes imposes gradients inside the cell [3] and, as we discuss below, may produce the polarization of the cell.
The experimental study of the interactions among MARCKS proteins, phospholipids and PKC has followed three different strategies: in vitro experiments wherein MARCKS interacts with vesicles composed by neutral and acidic phospholipids [12][13][14][15][16][17], in vitro experiments wherein MARCKS and PKC interacts with a Langmuir monolayer formed by a mixture of acidic and neutral phospholipids [18][19][20][21], and in vivo experiments with living cells subject to external perturbations [22][23][24][25]. In vitro experiments allow to control the phospholipid distributions on vesicles and in Langmuir monolayers and the concentrations of the proteins. They permit the characterization of the protein response to different concentrations of acidic phospholipids and the realization of dynamical experiments by the introduction of activated PKC [15,21]. A recent study for the interaction of MARCKS and PKC with a Langmuir monolayer [21] shows oscillations based on the nonlinear nature of protein-phopholipid interactions. In vivo experiments are more realistic, however, additional components in the cell may change the dynamics. The external http://www.epjnonlinearbiomedphys.com/content/2/1/1 introduction of a chemical signal produces variations of the Ca 2+ concentration, and can activate PKC. Experiments with living cells show the co-distribution of domains of high MARCKS concentration at the membrane with domains of acidic lipids [26] or domains of calmodulin [27]. The interaction of lipids and MARCKS proteins can lead to formation of acidic lipids covered by electrostatically bound proteins [28].
Experiments with vesicles as well as in-vivo experiments have motivated the derivation of different reaction-diffusion models of the interaction of MARCKS with membranes [29,30]. These models conserve the total number of MARCKS proteins, i.e., neither synthesis nor degradation of proteins are considered, because the time scales studied here is small in comparison with those of protein synthesis. Such reaction-diffusion models with mass conservation have been previously employed in the study of other systems, for example to explain how cell polarity emerges [31][32][33][34] to describe the protein system that controls cellular division in the bacterium E. Coli [35][36][37].
In this paper, we add the dynamics of calcium and PKC to a mass-conserved reactiondiffusion model for binding, phosphorylation and dephosphorylation of MARCKS proteins in living cells. The enzyme PKC is activated by Ca 2+ . The activated PKC binds at the membrane, where phosphorylates MARCKS proteins. In contrast with previous models [29,30], here, the parameter values were adapted to quantitatively reproduce the experimental results by Mogami et al. [23]. We compare the outcome of experiments in living cells after an increase of Ca 2+ concentration with the numerical results obtained with our model following an equivalent increase of calcium. Then, we perform numerical simulations of the spatio-temporal dynamics keeping a constant Ca 2+ concentration and the parameter values obtained by the comparison with experiments. Such calculations show the spontaneous formation of domains with high concentration of MARCKS proteins at the membrane. It suggests that the formation of domains of MARCKS proteins are likely to occur in living cells.
The manuscript is organized as follows. In the next section we describe the model for the PKC, MARCKS and Ca 2+ spatio-temporal dynamics and the numerical methods employed for the integration of the resulting model. Then, we obtain the parameters of the model from the literature and by comparing the resulting numerical simulations with the experimental results obtained from [23]. Once the parameter values are chosen, extended numerical simulations show the formation of domains of MARCKS proteins at the cellular membranes by two mechanisms. We analyze both mechanisms and calculate the phase diagram. Finally, we discuss the main results of the manuscript and, in particular, their biological implication.

Methods
We employ a reaction-diffusion model for the concentrations of MARCKS and PKC in the cell. First only the temporal dynamics of the concentrations described by ordinary differential equations is considered, then the complete model is derived taking into account transport processes and the spatial aspects of the dynamics.

Temporal reaction dynamics
The binding and unbinding processes of MARCKS and PKC proteins, together with the phosphorylation and dephosphorylation reactions of MARCKS proteins are studied. http://www.epjnonlinearbiomedphys.com/content/2/1/1

MARCKS dynamics
Phosphorylation and dephosphorylation of proteins are usually modeled by the Goldbeter-Koshland mechanism [10]. Using a quasi-steady approximation of the intermediate complexes, the original dynamics can be simplified to a couple of equations for the free and phosphorylated MARCKS (respectively [ MRK] f , and [ MRK] p ) using Michaelis-Menten rates [11]: where the total concentration of MARCKS is preserved. The terms k mrk b and k mrk u are the binding and unbinding rates of MARCKS. The rate of phosphorylation at the membrane and dephosphorylation in the cytoplasm are determined by the parameters k mrk pkc and k mrk ppa . The phosphorylation depends on the concentration of PKC at the membrane [ PKC] m and, therefore, indirectly on Ca 2+ . The dephosphorylation is assumed to be independent of [Ca 2+ ] [38] because MARCKS proteins can be dephosphorylated by the catalytic subunit of protein phosphatase 2A (PP2Ac), and protein phosphatase 1 (PP1c), which activation is independent of Ca 2+ [39]. A possible dependence of the phosphatase on Ca 2+ would induce a reciprocal regulation of the enzymes producing more complex dynamics [40].
The binding of MARCKS depends on the surface term:

PKC dynamics
We consider a simple binding and unbinding dynamics of the PKC enzymes at the membrane. The total concentration of enzymes is conserved giving rise to the next two coupled equations: [ PKC] c +k

Ca 2+ dynamics
Usually, [Ca 2+ ] is constant and stays at rest, however external signals or internal processes produce random spikes of large [Ca 2+ ]. In such case, the global dynamics of [ Ca 2+ ] is simply modeled by a relaxation dynamics to the rest state and random spike generation: where τ is the characteristic relaxation time to the rest state [ Ca 2+ ] o . This time is controlled by the concentration of buffers in the cell. The term I(t) represents the random spikes. To incorporte the spikes we increase instanstaneously the calcium concentration by [ Ca 2+ ] ∼ 0.03 − 0.1 μM. http://www.epjnonlinearbiomedphys.com/content/2/1/1

Reaction-diffusion equations
MARCKS and PKC are hydrophilic proteins and they transport to the membrane by diffusion. Here, for simplicity, we neglect effects of molecular crowding [41] and To account for this we divide the system in two regions: one corresponding to the membrane and its immediate vicinity, where all the translocation processes between membrane and cytoplasm take place, and the other to the bulk of the cytoplasm. For more details see the subsection on numerical methods below. The set of reaction-diffusion equations for the concentrations at the translocation zone reads: The equations describing the dynamics in the bulk of the cytosol read: where the only chemical transformation is the dephosphorylation of MARCKS. We restrict our study to global signals of [ Ca 2+ ] given by Eq. (4). A full description of the spatio-temporal dynamics of [ Ca 2+ ] requires request a detailed description of the mechanisms of the intracellular [ Ca 2+ ] dynamics, which is beyond the purpose and scope of this manuscript.

Parameter values
Some of the parameters of the model can be obtained or estimated from the literature. Typical values of the diffusion coefficient of proteins in cells are around D ≈ 1 − 10 μm 2 s −1 [42,43]. The diffusion coefficient of a kinase has been explicitly measured [44], and we use this typical result as characteristic diffusion coefficient of PKC D pkc c ≈ 5 μm 2 s −1 . http://www.epjnonlinearbiomedphys.com/content/2/1/1 Although the molecular weight of MARCKS (32Kda) is smaller than the weight of PKC (80Kda), its apparent weight is similar (80Kda) [45]. Therefore we employ the same value for the diffusion coefficient of MARCKS D c ≈ 5 μm 2 s −1 .
Membrane diffusion of proteins is smaller than the cytoplasmatic diffusion [46]. In particular, for mammalian cells, membrane diffusion coefficient of proteins has been estimated to be D m ≈ 0.1 μm 2 s −1 [47].
Calcium rest concentration has been previously estimated [48]. Complete dephosphorylation of MARCKS protein is obtained after less than 3 minutes in experiments with neutrophils [9], therefore, we infer a characteristic time of around 20 seconds. The total concentration of MARCKS [ MRK] = 10μM [4] is obtained from experiments with living cells. The coverage concentration of MARCKS is estimated for a typical cell in [29,30] and is of the same order than the total concentration.
Following [49,50] we take [ PKC] 0 ≈ 0.5 μM as typical value of PKC concentration. However, some groups have considered lower [51] or higher [52] concentration of PKC, in the range of The rest of model parameters are estimated by the fitting of experimental curves shown in [23]. The list of parameter values is shown in Table 1.

Numerical methods
For a better comparison with the experiments shown in [23], we consider a circular domain representing a cell with radius R = 10 μm. Note that the use of threedimensional environments with two-dimensional membranes does not qualitatively change the results, for more details see [29]. In the center we introduce a circular passive nucleus which does not interact with the dynamics, thus we consider non-flux boundary conditions in the internal radius (R i = 2.5 μm) of the integration domain.
We discretize the resulting corona using polar coordinates. We employ finite differences, adapted to polar coordinates, for the calculation of the diffusion terms. While the radial discretization is constant and corresponds to r = 0.125 μm, the angular spatial discretization depends on the radial coordinate.
The exterior layer, corresponding to the external radius of the cell, corresponds to the thin layer where the binding and unbinding dynamics to the membrane during time t occurs. The discretization of the system implies that this external layer is 0.125 μm thick, much larger than the real thickness of the cellular membrane, however, the use of finite differences permits the definition of this interaction length, for more details see [29]. Within this discretization, we employ the membrane concentrations [ MRK] m and [ PKC] m as volume concentrations.
We employ an explicit Euler method with t = 0.2 ms for the numerical integration of the equations.

Results
First we employ experimental data obtained from [23] to build a model by fitting the parameter values. Then, we study by numerical simulations of the resulting model different types of symmetry breaking mechanisms.

Building a quantitative model
We consider experiments on MARCKS and PKC interaction in INS-1 cells, an insulinsecreting cell [23] to obtain the values of the parameters.

Parametrizing PKC and MARCKS binding
The enzyme PKC is activated by Ca 2+ [38]. Once the enzyme is activated it can bind to the membrane. The PKC dependence on Ca 2+ seems to follow a Michaelis-Menten dynamics [53,54], however, the explicit form is derived from the experiments on living cells [23]. In such experiments there is a sharp transition between the concentrations of low active PKC and high active PKC at [ Ca 2+ ] = 400 nM. Low Ca 2+ concentrations ([ Ca 2+ ] < 400 nM) do not activate PKC and high Ca 2+ concentrations ([ Ca 2+ ] > 400nM) activate PKC and produce the binding of PKC at the membrane [23]. The relaxation times (τ ) observed in experiments are distributed between 5 and 15 seconds. For the numerical simulations shown here we use the value τ = 8 s. A close estimation is obtained in a similar context in [55].
In Figure 2 the responses of both, PKC and MARCKS, after a spike of Ca 2+ are shown simultaneously. Experimental data has been extracted form [23]. The ratio between the concentrations after and before the Ca 2+ spike is plotted for the experimental and numerical cases. A spike of [Ca 2+ ] decreases the concentration of PKC in the cytoplasm, see  Table 2.
The binding of the PKC enzymes at the membranes causes the unbinding of MARCKS proteins. The processes of binding and phosphorylation of MARCKS are fitted from the experiments to give rise to the similar behavior shown in Figure 2. We employ the comparison between these curves to fit the parameters for the binding (k mrk b ) and unbinding (k mrk u ) of MARCKS, and the phosphorylation rate (k mrk pkc ).

Dynamics of MARCKS and PKC
The absolute responses of MARCKS and PKC proteins to the Ca 2+ spike are shown in Figure 3. The concentration of PKC at the membrane increases, see Figure 3b, and phosphorylates some membrane MARCKS proteins. It produces an increase of phosphorylated MARCKS, see Figure 3c, and the subsequent dephosphorylation in the cytoplasm. Note that the low dephosphorylation rate produces a delay in the increase of free MAR-CKS proteins on the cytoplasm with respect the phosphorylated MARCKS. It increases the recovery time of MARCKS to the initial state with respect to the PKC dynamics. Such increase is shown in Figure 2 and has been observed in different experimental realizations [23].

Response to a set of [Ca 2+ ] spikes
The model can reproduce other experimental results. In Figure 4 we show the response of a different cell of the same type to a set of Ca 2+ spikes. Keeping the same parameter values obtained from Figure 2 the model reproduces the PKC and MARCKS responses, where we have only modulated the amplitude of the spikes of calcium. Whereas the response of MARCKS in the numerical simulation is comparable to the experiments, the comparison of the PKC response is not completely satisfactory. Tuning the value of some kinetic parameters or total concentrations may improve the results. However, we obtain a qualitative agreement without any further tunning.

Pattern formation in living cells
After the characterization of the model parameters from the above simulations, we perform large simulations for the whole cell, keeping the parameters constant and [Ca 2+ ] at rest.

Spontaneous domain formation
The evolution of the membrane concentration of MARCKS and PKC proteins are shown in Figure 5. The spatio-temporal plot, see Figure 5a domains at the membrane. First, the proteins are homogeneously distributed at the membrane, see Figure 5b. However, small spatial differences in the concentrations induced by small random perturbations lead to a long-wave instability. The spatial differences on the concentration increase with time and give rise to domains.
[PKC] accumulates at the membrane outside the domains of MARCKS proteins. Cells with larger number of transitory domains can be observed upon change of the values of some parameters or the cell size. These domains may be relatively stable, however, coarsening of domains is observed (results not shown) and, finally, a single large domain survives at the membrane.
To characterize the growing rate of the domains, the relative difference ( [ MRK] m ) between the local maximum and the minimum of [ MRK] m is plotted in Figure 5c. This difference corresponds to and saturates after an exponential growth. This growth defines a characteristic temporal scale for the growing rate (τ g = 285 s).

Bistability-induced domain formation
A second mechanism of pattern formation is associated with the bistability of two different homogeneous steady states. An inhomogeneous initial perturbation produces the simultaneous appearance of two stable solutions, which compete and finally develop a stationary profile.
A numerical simulation with higher MARCKS and PKC concentrations is shown in Figure 6a,b. The initial condition distributes the protein concentrations inhomogeneously along the cell giving rise to two different regions with high and low protein concentrations. It produces a polar orientation of the cell.
Small perturbations of the homogeneous initial state as in the previous section do not produce the formation of the domain, see Figure 6c,d. The distribution of MARCKS and PKC at the membrane arrives to two completely different states depending on the initial condition.

Change on MARCKS and PKC concentrations
In order to check the robustness of the pattern formation mechanisms, we have performed a systematic numerical study of the conditions of appearance of the membrane domains. We have chosen the total concentrations of MARCKS and PKC proteins as control parameters because they may be controlled in experiments [51]. The two types of pattern formation appearing in the model are shown in Figure 7. Simulations done with an small random perturbation on the initial condition can produce the formation of domains, light gray region in Figure 7. A single realization is shown in Figure 5. On the other hand, dark gray region corresponds to simulations which homogeneous solution is stable under small perturbations, but unstable if the perturbation is strong enough, see Figure 6.
Large concentrations of MARCKS saturate the membrane, above gray regions in Figure 7, whereas low concentrations do not allow the formation of domains, below gray regions in Figure 7. The absence of PKC precludes domain formation, in contrast, at fixed MARCKS concentration, large concentrations of PKC phosphorylate the membrane too fast and also suppress the domain formation.
The dashed line in Figure 7 corresponds to the cover concentration of MARCKS [ MRK] o . This quantity is a threshold for MARCKS concentration at the membrane.

Change on the diffusion coefficients
The two mechanisms responsible of domain formation are controlled by the values of the diffusion coefficients. In Figure 8 we study the dependence of the pattern formation http://www.epjnonlinearbiomedphys.com/content/2/1/1

Discussion
Formation of MARCKS domains has been observed in experiments in hippocampal neuron growth cones [26]. They are correlated with domains of acidic phospholipids as PIP 2 or PS. Experiments with vesicles show that phospholipid domains do not form in absence of MARCKS. The introduction of MARCKS proteins induces the formation of domains [17]. Our results are obtained without considering the dynamics of the phospholipids. The electrostatic interaction of the phospholipids with MARCKS may induce the formation of phospholipid domains following the protein domains. Domains of MARCKS proteins at the membrane of smooth muscle cells co-distributed with calmodulin have been observed in experiments [27]. We propose two different mechanisms of pattern formation for MARCKS protein at membranes: A long wave instability and a bistability-related mechanism. The long-wave instability is responsible for the pattern formation for the parameter values fitted from [23]. Increasing the concentration of MARCKS and PKC produces the appearance of domains induced by the second mechanism.
The formation of domains is robust under change of parameter values. Figure 7 shows the window of MARCKS and PKC concentrations which permit pattern formation. The relation among the values of the two diffusion coefficients (D m and D c ) which produce domain formation is shown in Figure 8. For a given value of D m pattern formation happens above a critical size of the cell.
The predicted domains may have a regulatory function by making the protection of the acidic phospholipids through membrane-bound proteins more efficient. In this way a local accumulation of proteins may prevent the effect of PKC inside the domain, which produces the phosphorylation and membrane-unbinding of MARCKS and the subsequent hydrolysis of the unprotected PIP 2 , which in turn induces a release of Ca 2+ [56].
The relation between MARCKS proteins and the induction of Ca 2+ spikes is beyond the scope of the present work. Here we have assumed that Ca 2+ follows a simple relaxation dynamics, see Eq. (4). However, it is known that Ca 2+ has a spatio-temporal dynamics, and traveling waves are generated by a calcium-induced calcium release mechanism [11]. The use of an adequate model for the Ca 2+ dynamics combined with the generation of Ca 2+ produced by the hydrolysis of the PIP 2 may provide a feedback loop between the processes at the membrane and in the interior of the cell, giving rise to a suitable mechanism of cell signaling control.
Possible stochastic effects due to a low number of proteins have been neglected here. The large concentration of MARCKS proteins in cells permits the use of deterministic dynamics. The concentration of PKC is smaller and stochastic effects may become relevant. The use of a stochastic model may enhance domain formation because the initially small perturbations necessary for pattern formation will naturally appear in a stochastic framework.
Living cells are three-dimensional. The model studied here is, however, a twodimensional approximation. The two-dimensional view is, however, sufficient to calculate proteins concentrations at the membrane. The realization of three-dimensional simulations would increase the computational time and the complexity of the numerical methods substantially.
The modelling results obtained in this work were compared to the dynamics of MAR-CKS and PKC in experiments with living insulin-secreting cells. However, methods and http://www.epjnonlinearbiomedphys.com/content/2/1/1 results are equally relevant for other systems involving membrane proteins. For example, the growth-associated protein 43 (GAP 43) and the cortical cytoskeleton associated protein (CAP 23) are proteins with similar properties as MARCKS [4]. They form also domains at the membrane of living cells [26] based on the interaction with acidic phospholipids and PKC [4]. On the other hand, the myelin basic protein (MBP), one of the most abundant proteins in the nervous system, has a structure analogous to MARCKS [57], electrostatically interacts with the PIP 2 of the membranes and responds to elevated Ca 2+ signals [58]. Furthermore, phosphorylation by PKC and other kinases regulate multiple processes in living cells, e.g. the formation of polarity of cells induced by PAR proteins [59], the cyclic dynamics of Rho GTPases [31] or the regulation of the cell division of E. Coli controlled by the Min proteins [37]. Hence, similar extensions of the Goldbeter-Koshland mechanism may be applicable to a large variety of biological systems.

Conclusions
In this paper, we have formulated a simple model for the binding and unbinding of MAR-CKS and PKC to the cellular membrane. Phosphorylation of MARCKS by PKC occurs at the membrane, while dephosphorylation of MARCKS in the cytosol. The parameter values of the model have been fitted from experiments shown in [23], wherein the binding dynamics of PKC and MARCKS are investigated. The resulting model with fitted parameters predicts the formation of protein domains at the membrane of the cell.
We employ an extended version of the classical Goldbeter-Koshland mechanism of protein phosphorylation and dephosphorylation. The separation of the two processes in two different compartments of the cell (phosphorylation at the membrane and dephosphorylation in the cytoplasm) combined with the difference of the diffusion coefficients at the membrane and in the cytoplasm (D m D c ) produces the formation of domains. Another important factor in the model is the conservation law included in the system.
In summary, we have developed a simple model of binding, phosphorylation and desphosphorylation for MARCKS protein. The main prediction of our numerical simulations is the spontaneous appearance of domains of high concentration of membrane proteins. The interaction of this protein domains with phospholipids at the membrane, in particular PIP 2 may have important consequences in the process of formation of IP 3 , and, therefore, in the control of the signaling mechanisms of Ca 2+ .