Modeling of regulatory loops controlling galactolipid biosynthesisin the inner envelope membrane of chloroplasts

Eric Maréchal, Olivier Bastien n

Laboratoire de Physiologie Cellulaire et Végétale, UMR 5168 CNRS-CEA-INRA-Univ., Grenoble Alpes, CEA Grenoble, 17 rue des Martyrs, 38054,Grenoble Cedex 09, France

H I G H L I G H T S

� We developed a mathematical model for plastid galactolipid synthesis by PAP and MGD1.� PAP is inhibited by diacylglycerol and MGD1 is activated by phosphatidic acid.� Phosphatidic acid is predicted to be channeled to phosphatidylglycerol synthesis.� An import of diacylglycerol is predicted to support galactolipid synthesis.� The model gives a possible explanation for the loss of the prokaryotic pathway.

a r t i c l e i n f o

Article history:Received 2 May 2014Received in revised form8 July 2014Accepted 10 July 2014Available online 18 July 2014

Keywords:Galactolipid biosynthetic pathwayDynamical systemMetabolic control analysisMetabolic processes

a b s t r a c t

In Angiosperms, the biosynthesis of galactolipids involves enzymes localized in the inner envelopemembrane (IEM) of chloroplasts, including a phosphatidic acid phosphatase (PAP), dephosphorylatingphosphatidic acid (PA) into diacylglycerol (DAG), and MGD1, transferring a galactose onto DAG thusgenerating monogalactosyldiacylglycerol (MGDG). It has been shown that PA and DAG could besynthesized in the plastid via the so-called ‘prokaryotic’ pathway or imported from the endoplasmicreticulum via the ‘eukaryotic’ pathway. In vitro studies support the existence of (1) a negative regulationof the plastid PAP by DAG and (2) an activation of MGD1 by PA. We developed a mathematical model ofthe IEM galactolipid biosynthesis pathway to understand the properties of the system ruled by thepresence of these two regulatory motifs. We demonstrated that the design of the system implies that PAshould accumulate to levels that are not observed experimentally, regardless of its prokaryotic oreukaryotic origin. PA should therefore be used for other syntheses, such as that of phosphatidylglycerol.Whereas a massive influx of eukaryotic PA appears unlikely, an influx of eukaryotic DAG in the IEM issupported by simulations. The model also implies that DAG cannot transiently accumulate and that PAmainly acts as a signal switching the whole system on. Eventually, this analysis highlights the fact thatthe PAP enzyme could easily become dispensable and that the design of the system, with the tworegulatory motifs, could precede the loss of the PAP gene or activity in this pathway, a phenomenon thatoccurred independently in most clades of Angiosperms.

& 2014 Elsevier Ltd. All rights reserved.

1. Introduction

The bulk of chloroplast membrane glycerolipids is composed ofnon-phosphorous galactoglycerolipids, i.e. up to 50% monogalactosyl-dyacylglycerol (MGDG) and 30% digalactosyldiacylglycerol (DGDG).Other chloroplast lipids consist of sulfoquinovosyldiaclglycerol

(SQDG) and a single phospholipid, phosphatidylglycerol (PG) (Fig. 1).This composition is reminiscent of that found in cyanobacteria, fromwhich chloroplasts have derived, and contrasts with other subcellularmembranes, which are phospholipid-rich. The biosynthesis of MGDG,DGDG, SQDG and PG occurs within the two-membrane system thatdelineates chloroplasts, the plastid envelope, using two pools ofdiacyl-precursors, either produced in situ within chloroplasts orimported from other cell compartments (for recent review, Boudièreet al., 2012; Boudière et al., 2014). The precise nature of the importeddiacyl-precursors is still not resolved and as a consequence, thecorresponding metabolic fluxes within the cell are completelyunknown.

Contents lists available at ScienceDirect

journal homepage: www.elsevier.com/locate/yjtbi

Journal of Theoretical Biology

http://dx.doi.org/10.1016/j.jtbi.2014.07.0130022-5193/& 2014 Elsevier Ltd. All rights reserved.

n Correspondence to: Laboratoire de Physiologie Cellulaire et Végétale, Institut deRecherches en Technologies et Sciences pour le Vivant (iRTSV), CEA-Grenoble,17 rue des Martyrs 38054 Grenoble Cedex 9, France. Tel.: þ33 4 38 78 38 55;fax: þ33 4 38 78 50 91.

E-mail address: olivier.bastien@cea.fr (O. Bastien).

Journal of Theoretical Biology 361 (2014) 1–13

Schematically, the structure of membrane glycerolipids consistsof a three-carbone glycerol backbone (conventionally numberedsn-1,-2 and -3), esterified to fatty acids at positions sn-1 and sn-2and harboring a polar head at position sn-3 (Fig. 1). Two maintypes of diacyl-glycerolipid structures occur in plant cells and havedistinct acyl signatures. On the one hand, the ‘prokaryotic’ struc-ture is synthesized primarily in the inner envelope membrane(IEM) and contains a 16-carbon fatty acid at the sn-2 position, likecyanobacterial lipids (Fig. 2, pathway shown in red). On the otherhand, the ‘eukaryotic’ type is synthesized in the endoplasmicreticulum (ER) and contains an 18-carbon fatty acid at the sn-2position (Fig. 2, pathway shown in blue). Some plants, such asArabidopsis thaliana, have both prokaryotic- and eukaryotic-typeMGDG; they are called “16:3 plants”. Other plants, such as Pisumsativum, have only eukaryotic-type MGDG and are called “18:3plants” (Heinz and Roughan, 1983). The production of chloroplastMGDG relies therefore on an important fraction of substrateoriginating (and imported) from the ER. Regardless of the presenceof MGDG of prokaryotic or eukaryotic structure, DGDG is mainly ofthe eukaryotic type (Fig. 1). Galactoglycerolipids are assembled inthe chloroplast envelope following the pathway shown in Fig. 2(for review, Block et al., 2007; Boudière et al., 2014).

MGDG is synthesized by MGDG synthases (also called MGDs),which transfer galactose from a uridine-diphospho-galactose(UDP-Gal) donor onto a DAG acceptor:

MGD1 activity in the inner envelope membrane : DAGþUDP�Gal-MGDGþUDP

In Arabidopsis, there are three MGDG synthases, i.e. MGD1,MGD2 and MGD3. MGD1 is a MGDG synthase of type A, harboringa cleavable chloroplast transit peptide and targeted to the IEM(Awai et al. 2001). MGD2 and MGD3 belong to type B, without anytransit peptide, and localized in the outer envelope membrane(OEM) (Awai et al., 2001). MGD1 is the major (if not unique)isoform acting for the massive production of galactolipids neces-sary for the development of photosynthetic membranes (Jarvis etal., 2000; Awai et al., 2001; Kobayashi et al., 2007). MGD2 andMGD3 are expressed in non photosynthetic tissues and are alsoinvolved in some specific physiological contexts, such as theresponse to phosphate shortage (for review, Benning and Ohta,2005; Boudière et al., 2012).

Here, we shall focus on MGDG synthesis catalyzed by MGD1, inthe IEM and in normal developmental and physiological condi-tions. MGD1 is able to utilize either prokaryotic DAG synthesizedde novo in the IEM or eukaryotic DAG derived from ER

phosphatidylcholine (PC) (Maréchal et al., 1994a, 1994b; Miègeet al., 1999; Block et al., 2007) (see position of MGD1 in Fig. 2).

Concerning the generation of prokaryotic precursors for MGD1,PA is first produced in the IEM by the stepwise acylation ofglycerol-3-phosphate (Frentzen et al., 1983). It is then depho-sphorylated by a PA phosphatase (PAP), generating the prokaryoticDAG substrate by the following reaction:

PAP activity in the inner envelope membrane : PA-DAGþphosphate

The activity of the PAP enzyme in the IEM has not beenanalyzed yet using the totally purified enzyme or a recombinantsystem. It has only been analyzed in native fractions of spinachchloroplast envelope, using a synthetic PA substrate with 18-carbon at both sn-1 and sn-2 positions (dioleoyl-PA), i.e. aeukaryotic substrate (Malherbe et al., 1992). The relative specificityof PAP for prokaryotic versus eukaryotic PA is not known. Thisenzyme can provide precursors for the synthesis of MGDG, DGDGand SQDG. Note that the proportion of SQDG in plastid membranesis far lower than that of galactolipids and that the genetic knockout of its synthesizing enzyme has no critical impact on Arabi-dopsis phenotype, except when plants grow in nutrient shortage(Yu et al., 2002). We therefore designed our model without theSQDG synthesizing branch. Interestingly, in spinach envelopemembranes, PAP was shown to be inhibited by DAG, its product(Malherbe et al., 1992), regulating its activity by a retro-controlloop. This remarkable feature has been implemented in the modelwe developed here.

Alternatively, PA can be converted into cytidine diphosphate-diacylglycerol (CDP-DAG), a critical precursor for the synthesis ofPG (Xu et al., 2006), which is mostly of prokaryotic type (Fig. 1).Genetic engineering experiments have shown that eukaryotic-type PA could be also used for the synthesis of plastid PG (Fritz etal., 2007). It is therefore difficult to understand why, if eukaryoticPA is massively imported into the IEM, it would not be used for PGbiosynthesis. Without any other evidence, these studies suggestthat some kind of channeling might occur to direct prokaryotic-PAspecies to PG.

Concerning the origin of the eukaryotic structure, previousstudies have shown that PC could provide a diacyl-substrate forMGDG synthesis after the action of either non-specific phospho-lipase C (PLC) enzymes (NPC4 and/or NPC5, Nakamura et al., 2005;Nakamura et al., 2009) or phospholipase D (PLD) enzymes (PLDζ1and/or PLDζ2, Cruz-Ramírez et al., 2006; Oropeza-Aburto et al.,2012) (see Fig. 2). On the one hand, the action of PLCs releases DAG

CH2R1COO

CH

CH2

R2COO

P

CH2R1COO

CH

CH2

R2COO

P

CH2R1COO

CH

CH2OH

R2COO

CH2R1COO

CH

CH2OH

R2COO

OCH2OH

CH2R1COO

CH

CH2O

R2COO

OCH2OH

OCH2OH

CH2R1COO

CH

CH2O

R2COO

OCH2O O

CH2OH

CH2R1COO

CH

CH2O

R2COO

OCH2O O

CH2OH

OCH2O O

CH2OH

CH2R1COO

CH

CH2O

R2COO

CH2R1COO

CH

CH2

R2COO

P

CH2OH

CHOH

H2C

CH2R1COO

CH

CH2

R2COO

P

CH2OH

CHOH

H2C

phosphatidic acid(PA)

diacylglycerol(DAG)

sn-1

sn-2

sn-3

monogalactosyldiacylglycerol(MGDG)

digalactosyldiacylglycerol(DGDG)

phosphatidylglycerol(PG)

Fig. 1. Chloroplast galactoglycerolipids (monogalactosyldiacylglycerol, MGDG and digalactosyldiacylglycerol, DGDG) and their biosynthetic precursor (diacylglycerol, DAG)and regulators (phosphatidic acid, PA and phosphatidylglycerol, PG). PA is also an upstream precursor for the biosynthesis of DAG and PG.

E. Maréchal, O. Bastien / Journal of Theoretical Biology 361 (2014) 1–132

with the appropriate eukaryotic structure. The subsequent transferof DAG to the plastid envelope is supported by in vitro experimentsusing isolated chloroplasts (Andersson et al., 2004), but it has notbeen demonstrated yet at the cellular level. On the other hand, theaction of PLDs releases PA, which is possibly imported to the IEM,where it could be converted into DAG by the plastid PAP describedabove (Fig. 2). PA import is supposed to occur mainly via an ATP-binding cassette complex, the TGD transporter, made of threesubunits, TGD1, TGD2 and TGD3 in the IEM, associated to a fourthcomponent, TGD4 in the OEM (Xu et al., 2010; Wang et al., 2012;Wang et al., 2013). It has also been proposed that phospholipidsmight be transferred to the chloroplast (Mongrand et al., 2000),prior their hydrolysis by phospholipases, releasing PA or DAG. PLCand PLD reactions are summarized below:

PLC activity in endomembranes or at the surface of

the chloroplast : PC-DAGþphospho�choline

PLD activity in endomembranes : PC-PAþcholine

The MGD1 enzyme from Arabidopsis converts prokaryotic- andeukaryotic-type DAG into MGDG with the same specificity, i.e.identical Michaelis/Menten constant (KM) and identical maximalvelocity (VM) (Awai et al., 2001), by contrast with the spinachenzyme which has its highest affinity, i.e. lowest KM, for eukaryoticDAG (Maréchal et al., 1994a). Interestingly, the MGD1 activity wasshown to be regulated by anionic lipids, PG and PA (Dubots et al.,2010) (Fig. 1). This feature has also been used in the proposed model.

Take together, in vitro studies support therefore the existence oftwo ‘regulatory motifs’ in the biosynthetic pathway leading toMGDG synthesis in the IEM, firstly a negative regulation of theplastid PAP by DAG (Malherbe et al., 1992) and secondly anactivation of MGD1 by PA (Dubots et al., 2010). Given the twopossible sources of DAG and PA, synthesized in the IEM or divertedfrom the ER, the tuning of MGDG synthesis by these precursorsappears as one of the most critical regulatory processes ofmembrane lipid homeostasis at the whole cell level. This regula-tory systemwould then allow the selective utilization of PA or DAGas a eukaryotic precursor for MGDG synthesis and enable thesophisticated adjustment of chloroplast galactolipids by inter-mediates reflecting the metabolic statuses in both the organelleand the cytosol.

Nevertheless, it is still unknown how the existence of the tworegulatory motifs revealed by experimental studies (PAP inhibitionby DAG and MGD1 activation by PA) could influence the equili-brium concentration of lipids in the IEM and subsequently in thethylakoids, and the magnitude of the produced MGDG flux.Concentration and flux modeling and dynamic simulations mighthelp resolving the presence of these regulatory motifs, butdifficulties are raised by the lack of data on the following meta-bolic parameters:

– the PA and the DAG fluxes into the IEM are unknown;– the PAP and the MGD1 concentration into the IEM are

unknown;– the concentration of PA in the IEM of Arabidopsis thaliana wild

type chloroplasts can be deduced from that of thylakoids about

PA DAG

MGDGMGD1

TGD

DGD1

PA

MGDGDAG

thylakoids

stroma

IEM

OEM

ER +endomembranes

cytosol

chloroplastenvelope

ATS1, ATS2 DGDG

DAGPC

PLD PLC

?

?

DGDG

?

??

MGDG, DGDG

PAP

Fig. 2. Simplified scheme of the transport and synthesis of galactoglycerolipids in a chloroplast. In this representation, only MGD1 and DGD1, the main isoforms of MGDGsynthase and DGDG synthase in green tissues are shown. MGD1 is localized in the inner envelope membrane (IEM), whereas DGD1 is in the outer envelope membrane(OEM). MGD1 transfers a first galactose from a UDP-Gal donor onto a diacylglycerol (DAG) acceptor. DGD1 transfers a second galactose from UDP-Gal onto MGDG. Twomolecular structures can be synthesized, either a “prokaryotic” one, following a pathway located in the chloroplast (red arrows) or a “eukaryotic” one, following the importof precursors from the endoplasmic reticulum (blue arrows). The synthesis of “prokaryotic” MGDG is based on the upstream synthesis of phosphatidic acid (PA), by thestepwise action of glycerol-3-phosphate sn-1 acyltransferase (ATS1) and lyso-PA sn-2 acyltransferase (ATS2) in the IEM. PA is then dephosphorylated by a PA phosphatase(PAP) of the IEM. Alternatively “eukaryotic” PA can be produced in the endoplasmic reticulum (ER) or other compartments of the endomembrane system, by the hydrolysis ofphosphatidylcholine (PC) by phospholipase D (PLD) enzymes, such as PLDζ1 or PLDζ2. PA has an affinity for subunits of the TGD transporter and it has been proposed that itcould be imported to the IEM by this multiprotein complex. The TGD transporter is an ATP-binding cassette complex, made of three subunits, TGD1, TGD2 and TGD3 in theIEM, associated to a fourth component, TGD4 in the OEM (Xu et al., 2010; Wang et al., 2012; Wang et al., 2013). DAG with a “eukaryotic” structure can also be produced inendomembranes by the hydrolysis of PC by phospholipase C (PLC) enzymes, such as NPC4 or NPC5. DAG might then be imported to the IEM via an unknown transporter,although this import has not been unambiguously demonstrated. In non-green tissues or when plants grow under a shortage of phosphate, MGD2, MGD3 and DGD2 areexpressed and locate in the OEM: these enzymes, which are not primarily involved in the massive production of lipids for the membrane expansion of thylakoids are notconsidered in the present study. Unknown trafficking systems are indicated by question marks (?). (For interpretation of the references to color in this figure legend, thereader is referred to the web version of this article.)

E. Maréchal, O. Bastien / Journal of Theoretical Biology 361 (2014) 1–13 3

0.08 mol% (Dubots et al., 2010; Fritz et al., 2007), but theconcentration of DAG is unknown. The experimental processto obtain purified chloroplasts and then purify envelope mem-branes, triggers the activation of a galactolipid:galactolipidgalactosyltransferase (GGGT, Dorne et al., 1982 – now knownto correspond to SFR2, Moellering et al., 2010), that generatesartifactual amounts of DAG in envelope membranes. An experi-mental procedure has been developed to inactivate the GGGTenzyme by treatment with a protease (Dorne et al., 1982), andbased on the numerous preparations of IEM fractions in ourlaboratory, DAG is a very minor component, far below 1 mol%,used here as an upper limit. Besides this imprecise upperbound, no measure of DAG in chloroplast membranes iscurrently considered as accurate.

In addition to these unknown parameters, the pathway we aimto model here is not a classical series of enzymatic reactions in awater solution, like the majority of metabolic pathways occurringwithin the cell. PAP and MGD1 are embedded in a biologicalmembrane, one might consider as a planar surface, and theirlipidic substrates and products are constrained to move laterally,or transversally, within this membrane environment.

The model we develop here takes into account all thesedifficulties and represents, to our knowledge, one of the firstsimulation of a membrane lipid biosynthetic module.

2. Methods – design of the model

2.1. The kinetic models

2.1.1. The PAP kinetic modelThe IEM phosphatidate phosphatase catalyzes the dephosphor-

ylation of PA into DAG that remains in the IEM and releases awater-soluble phosphate product, which does not inhibit theenzyme activity at physiological concentrations and is thereforenot considered in our model. Using isolated envelope membranes,Malherbe et al. (1992) demonstrated that the PAP activity wasreduced, when the membrane was enriched in DAG. More pre-cisely, they demonstrated that DAG was a powerful competitiveinhibitor of the reaction. The conclusion was also extended toisolated intact spinach chloroplasts. As a consequence, the in situPAP activity can be modulated by its product, via a classical retro-control loop. Starting from this conclusion and after havingconsidered the Lineweaver–Burk plots presented in this study(Malherbe et al., 1992), we can derive a model for PAP kinetics,following Cornish-Bowden (2004):

υnPAP ¼VM;PAP ½PA�

KM;PAPð1þð½DAG�=KI;PAPÞÞþ½PA� ð1Þ

where υnPAP, VM,PAP, KM,PAP, KI,PAP, [PA] and [DAG] are the enzyme'svelocity, the maximum enzyme velocity, the Michaelis–Mentenconstant (Michaelis and Menten, 1913), the DAG-inhibition con-stant, and the PA and DAG concentrations, respectively. Theenzyme's velocity υnPAP is usually given as a ‘specific activity’, i.e.in moles of lipids min�1 mg of protein�1, allowing an estimate inmolecules of lipids s�1 molecule of enzyme�1 (Malherbe et al., 1992,Dubots et al., 2010). When a substrate or inhibitor/activator is anamphipatic compound embedded within the membrane, its sur-face concentrations is given in mole fraction, calculated by theamount of this compound (in moles) divided by the total amountof all constituents in the membrane. This value is dimensionless(Cohen, 1996) and can be also expressed in mol% (¼ molefraction�100). While this convention is experimentally justifiedfollowing the surface dilution model (Dennis et al., 1981; Maréchal

et al., 1994a; Salinas et al., 2011), it implies both dramatictheoretical problems (units of the maximum velocity and of theMichaelis–Menten constant, apparent or not, must be redefined inorder to have a kinetic equation with consistent dimensions in theleft and right members) and numerical problems (coupling equa-tions of this form will lead to units’ inconsistencies).

Based on the studies by Malherbe et al. (1992), we can extractnumerical values for VM,PAP, KM,PAP and KI,PAP parameters (Table 1).A factor called W was introduced to correct the dimension of theenzymatic activity so that it fits with the lipid concentration pertime unit, W¼number of MGD1 molecules/total amount of lipidmolecules (the choice of the number of MGD1 molecules isjustified by the non-dimensionalization procedure, see below).

υPAP ¼WVM;PAP ½PA�

KM;PAPð1þð½DAG�=KI;PAPÞÞþ½PA� ð2Þ

where υPAP is a corrected expression of the enzyme's velocity,expressed in time unit�1.

2.1.1.1. The MGDG synthase 1 kinetic model. The IEM monogalacto-syldyacylglycerol synthase 1 catalyzes the transfer of a galactosylgroup from a water-soluble UDP-Gal molecule onto DAG,generating MGDG that remains in the IEM and UDP, a water-soluble product. UDP-Gal and UDP are not considered in thepresent model, based on the assumptions that UDP-Gal is notlimiting and UDP, which has been shown to inhibit MGD1(Maréchal et al., 1994a), is not in sufficiently high concentrationto exert a significant effect on the system.

It has been known for long that MGD1 is activated in presenceof anionic lipids (Covés et al., 1988) but the precise analyses of PAand PG activation have only been dissected recently (Dubots et al.,2010). MDG1 requires at least one of these two lipids to show anactivity, the most efficient activator being PA. One of the majorresults by Dubots et al. (2010), which is critical for the presentmodel, is the fact that “PA switches MGD1 on” and has a majoreffect on the maximum velocity measured for this enzyme, VM,

MGD1, but no detectable effect on the Michaelis–Meten constant,KM,MGD1. As a conclusion, MGD1 activation does not follow theclassical mechanism of an essential activation (which leads to themodulation of the KM with no VM variation) or mixed activation(which leads to KM and VM modulations) (Cornish-Bowden, 2004).In the same way, a hyperbolic activation should not be con-sidered as no activity could be observed in the absence of PA(Dubots et al., 2010). Taken together, experimental results supportthe binding of the activator to the substrate–enzyme complex.

Table 1Kinetic parameters of the developed model and numerical estimates based onin vitro experimental data. All parameters are defined in the text.

Parameters Estimated numerical values References

VM,PAP 4 Molecules of PA s�1.molecule of MGD1�1

Malherbe et al. (1992)

KM,PAP 6.0 mol% Malherbe et al. (1992)KI,PAP 0.7 mol% Malherbe et al. (1992)VM,MGD1 20 Molecules of DAG s�1.

molecule of MGD1�1Nishiyama et al. (2003)

KM,MGD1 1.0–5.0 mol% Maréchal et al. (1994)Awai et al. (2001)

KA,MGD1 1.5 mol% Dubots et al. (2010)h 1 Maréchal et al. (1994a)W Unknownα Unknownβ Unknownη Unknown

E. Maréchal, O. Bastien / Journal of Theoretical Biology 361 (2014) 1–134

We can construct the following model for MGD1:

υMGD1 ¼WVM;MGD1½DAG�

KM;MGD1þð1þðKA;MGD1=½PA�ÞÞ½DAG�ð3Þ

where υMGD1, VM,MGD1, KM,MGD1, KA,MGD1, [PA] and [DAG] are theenzyme's velocity, the maximum enzyme velocity, the Michaelis–Menten constant, the PA-activation constant, and the PA and DAGconcentrations, respectively. In this model, UDP-Gal concentrationis assumed to be non-limiting and constant and VM,MGD1 and KM,

MGD1 apparent Michaelis–Menten parameters (Cornish-Bowden,2004). From the enzymatic studies performed by Maréchal et al.(1994a), Nishiyama et al. (2003) and Dubots et al. (2010) we canextract numerical values for MGD1 parameters (Table 1).

2.2. Model of the chloroplast IEM galactolipid biosynthetic pathway

2.2.1. Design of the biological modelFrom the previous considerations, we constructed a model for

the galactolipid biosynthetic pathway occurring in the IEM, basedon an ordinary differential equation system. In this system,illustrated in Fig. 3, the α parameter represents the net influx ofPA, combining the in situ synthesis by the prokaryotic pathwayand the import of eukaryotic PA from the ER (see Fig. 2). Inaddition to DAG synthesis via the action of PAP, a parameter, η, wasintroduced to represent a putative influx of DAG from the OEM.Both α and η are assumed to be constant during the experimentaltime laps. Downstream, β([MGDG]) represents the exit flux ofMGDG from the system, exported to other membrane compart-ments, i.e. the OEM or the thylakoids (see Fig. 2). In the OEM,MGDG can serve as a substrate for the DGDG synthase activity.

d½PA�dt ¼ α� WVM;PAP ½PA�

KM;PAP ð1þð½DAG�=KI;PAP ÞÞþ ½PA�d½DAG�

dt ¼ WVM;PAP ½PA�KM;PAP ð1þð½DAG�KI;PAP ÞÞþ ½PA��

WVM;MGD1 ½DAG�KM;MGD1 þð1þðKA;MGD1=½PA�ÞÞ½DAG�þη

d½MGDG�dt ¼ WVM;MGD1 ½DAG�

KM;MGD1 þð1þðKA;MGD1=½PA�ÞÞ½DAG��βð½MGDG�Þ

8>>>><>>>>:

ð4Þ

The last equation in the model (4) is uncoupled from the firsttwo equations: it simply gives the relation between the concen-tration of MGDG, [MGDG], in mole fraction, with that of PA, [PA]and DAG, [DAG], once those have been determined. So we onlyneed to analyze the first two equations of this system.

A first observation is that the stationary “S” point for [MGDG],if there is one, is the value named [MGDG]S given by solvingαþη�β([MGDG])¼0. In general, β([MGDG]) is an increasing

function of [MGDG], being null when [MGDG]¼0. Three majorcases must be distinguished:

1. If β([MGDG]) is unbounded, then at least one [MGDG]S valueexists, for which αþη–β([MGDG])¼0.

2. If β([MGDG]) is bounded, then αþη–β([MGDG])¼0 will not besatisfied for large values of α and/or η.

3. If β([MGDG]) is not monotonic, then multiple stationary pointscould be allowed.

2.2.2. Expression in non-dimensional termsBefore analyzing the model (4), a necessary step is to express it

in non-dimensional terms (Murray, 2003; Segel, 1972). Two mainsadvantages come with non-dimensionalization: on the one hand(i) units vanish and so, the adjectives “small” and “large” have adefinite relative meaning, and on the other hand (ii) the number ofrelevant parameters is reduced to a few dimensionless groupings,which then determine the dynamics of the model (Murray, 2003;Segel, 1972). We thus introduced the following non-dimensionalquantities:

u¼ ½PA�KM;PAP

; v¼ ½DAG�KM;PAP

; τ¼ tWVM;PAP

KM;PAP; f ¼ α

VM;PAPW; d¼ η

VM;PAPW;

b¼ KM;PAP

KI;PAP; c¼ KM;PAP

KA;MGD1; g ¼ VM;MGD1

VM;PAP; r¼ KM;MGD1

KM;PAP:

Using u, v, τ, f, d, b, c, g and r, the equations expressed in(4) become

dudτ

¼ f � uð1þbvÞþu

dvdτ

¼ uð1þbvÞþu

� gvrþð1þð1=cuÞÞv þd

8>>><>>>:

ð5Þ

whereas the model was initially based on 9 parameters and3 variables, the non-dimensional model has only 6 parameters,which are pure numbers and 3 variables, u, v and τ. For example,if u⪡1, it simply means that the PA mole fraction is very smallcompared to the Michaelis–Menten constant of the PAP, alsoexpressed in mole fraction.

2.3. Metabolic control coefficients

We then analyzed the metabolic model shown in Fig. 2 usingthe conventional Metabolic Control Analysis (MCA) conceptsdeveloped by Kacser and Burns (1973), and Heinrich andRapoport (1974), based on previous works by Higgins (1963).Although it is rarely made explicit, MCA is a first order sensitivityanalysis in the vicinity of a structurally stable fixed point. MCAallows quantifying how metabolic variables, such as fluxes andspecies concentrations, depend on network parameters (Fell,2005). In the MCA framework, one is able to describe how networkdependent properties, called control coefficients, depend on localproperties, called elasticity coefficients (Heinrich and Schuster,1996). For a given metabolic network, e.g. metabolites Al, Am, An,…, an elasticity coefficient ευiAl

, is defined as the relative change of areaction rate of the i-th reaction (υi) in response to a relativechange of the l-th network species Al, that is to say

ευiAl¼ ∂ log υi∂ log Al

The two main control coefficients are on the one hand the flux andon the other hand concentration control coefficients. They mea-sure the relative steady state changes of a pathway flux (Φ) ormetabolite concentrations (Al, Am, An,…), in response to a relative

PA DAG MGDG

(-)

PAP MGD1

(+)

Fig. 3. Model of the regulatory motifs controlling MGDG biosynthesis by MGD1.In this simple model, PA can be produced in the IEM (prokaryotic PA) or importedfrom the ER (eukaryotic PA). At steady state, the influx of both prokaryotic andeukaryotic PA is expressed as a constant, α. DAG is produced in the IEM by theaction of the plastid PAP, or imported from the OEM or other subcellularcompartments. At steady state, the influx of eukaryotic DAG is expressed as aconstant, η. Once produced by MGD1, MGDG can leave the system, being eitherexported to the OEM to serve as substrate for the synthesis of DGDG, or desaturatedby FAD5 and other desaturases, and transferred to thylakoids. At steady, the flux ofMGDG generated by the system is expressed as a constant, β. Inhibition of PAPactivity by DAG and activation of MGD1 by PA are shown by arrows, (�) and (þ)respectively.

E. Maréchal, O. Bastien / Journal of Theoretical Biology 361 (2014) 1–13 5

change in a parameter, like an enzyme activity at the level of the i-th reaction (υi).

Flux control coefficients are defined by

CΦυi¼ ∂ log Φ

∂ log υi

Concentration control coefficients for the l-th species are definedby

CAlυi¼ ∂ log Al

∂ log υi

The flux control summation theorem (Kacser and Burns, 1973;Heinrich and Rapoport, 1974) states that metabolic fluxes aresystemic properties and that their control is shared by all reactionsin the system. In other words, if a single reaction changes itscontrol of the flux, this change is compensated by changes in thecontrol of the same flux by all other reactions (Cornish-Bowdenand Cárdenas, 1990). Given all the enzymatic reactions of thesystem (all possible values of enzymatic reactions i), the expres-sions for flux and l-th concentration control summation theoremsare ∑iC

Φυi¼ 1 and ∑iC

Alυi¼ 0. Relationships between elasticity

coefficients (ε) and control coefficients (C) are given by theconnectivity theorems. They highlight the close relationshipbetween the kinetic properties of individual reactions and thesystem properties of a pathway. Their expressions for flux andconcentrations are {∑iC

ΦυiευiAl

¼ 0} and {∑iCAnυiευiAm

¼ �1 if n¼m or∑iC

AnυiευiAm

¼ 0 if nam}, respectively. For linear pathway, connectiv-ity and summation theorems allow the complete determination offlux and concentration control coefficients since elasticities areknown. Nevertheless, once new features are introduced in thepathway, such as branch points, the number of enzymes exceedsthe number of internal metabolites and additional relationshipsbetween control coefficients must be found: these are the branchpoint theorems (Kacser, 1983; Fell and Sauro, 1985; Sauro et al.,1987).

2.4. Optimization procedures

Given f and d, the computation of the steady states for Eq. (5),both analytically and numerically, was performed using thesimplex algorithm (Dantzig, 1998) implemented in Matlab (TheMathWorks, 2013). For all f and d in a given interval (see Section 3.3.),steady states for non-dimension [PA] and [DAG] were computed andelasticity coefficients, fluxes and species control coefficients werededuced from these values (see Section 3).

3. Results

3.1. Stationary points

3.1.1. Existence of a stationary pointThe pathway shown in Fig. 2 has been modeled based on the

simple assumptions and concepts described above, eventuallyexpressed in dynamic terms in the system of Eq. (5). In this study,we only consider physiological situations at steady states, makingsense in homeostatic terms. If present, a steady state (us,vs) for(5) is solution of

Γðu; vÞ ¼ dudτ

¼ f � uð1þbvÞþu

¼ 0

Δðu; vÞ ¼ dvdτ

¼ uð1þbvÞþu

� gvrþð1þð1=cuÞÞv þd¼ 0

8>>><>>>:

ð6Þ

Considering that Γ(u,v) and Δ(u,v) are linearly independent(i.e. the Wronskian determinant of these two differential functions

does not vanish identically; Peano, 1889), the system (6) isequivalent to the following one, which is easier to express interms of equations:

Γðu; vÞ ¼ 0Γðu; vÞþΔðu; vÞ ¼Πðu; vÞ ¼ 0

(

Thus, the system that has to be solved is

Γðu; vÞ ¼ f � uð1þbvÞþu

¼ 0

Πðu; vÞ ¼ f � gvrþð1þð1=cuÞÞv þd¼ 0

8>><>>: ð7Þ

The Eq. Γ(u,v)¼0 defines a curve, v1(u), in the (u,v) plane withequation v1(u)¼(1/b)((1/f�1)u)�1/b (Fig. 4). A necessary condi-tion for the existence of a stationary point is given by the obviouscondition vsZ0, which required fo1. In dimensional terms, this isequivalent to the obvious condition: αoW.VM,PAP (Table 2). Whenthis condition applies and if stationary points exists, we havetherefore us41/(1� f).

The Eq. Π(u,v)¼0 (Fig. 4) defines a curve, v2(u), in the (u,v)plane with equation v2ðuÞ ¼ ðru=�ð1=cÞþððg=f þdÞ�1ÞuÞ. Again, anecessary condition for the existence of a stationary point is givenby the obvious condition vsZ0, which requires g4 fþd. Indimensional terms, this implies that αþηoW.VM,MGD1 (Table 2).Combining the inequality g4 fþd with the previous fo1, we thendeduce the following necessary condition for a steady state:fþdomin(dþ1, g).

As v2(u) is strictly decreasing on �ð1=cððg=f þdÞ�1ÞÞ; þ1½ andv1(u) is strictly increasing on the same interval, they intersect onlyin one point and so, there is only one steady state given bythe system (6) (intersection of curves shown in Fig. 4). Thestationary point is explicitly given by us ¼ ð�ðbrcþX�YcÞ�ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiðbrcþX�YcÞ2þ4XYc=2XYc

qÞ and vs ¼ ð1=bÞððð1=f Þ�1Þus�1Þ,

where X ¼ ð1=f Þ�1 and Y ¼ 1�ðg=dþð1=Xþ1ÞÞ (see Appendix A).In biological terms, given α a value of an influx of PA, and η an

influx of DAG, in the pathway shown in Fig. 2, only one steady stateexists and it is given by the above formula. We can thereforesimulate, compare and analyze various reasonable scenarios withvariable influxes of PA and DAG and determine whether thesescenarios are plausible, when compared to natural contexts.

0 0.5 1 1.5 2 2.5 30

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

0.18

0.2

Γ(u,v)=0Δ(u,v)=0Δ'(u,v)=0

u

v

Γ(u,v)=0Δ(u,v)=0Π(u,v)=0

Γ(u,v)=0Δ (u,v)=0Π(u,v)=0

Fig. 4. Existence of stationary states for the model. The existence of steady state(s)(us,vs) in system (5) is solution of Γ(u,v), Δ(u,v) and Π(u,v) being null (see text).Based on the simulation of Γ(u,v), Δ(u,v) and Π(u,v) null clines, the intersection of Γ(u,v)¼0 with Δ(u,v)¼0 (or Γ(u,v)¼0 with Π(u,v)¼0) allows the identification ofonly one steady state. Shown simulation was performed with W¼0.2, α¼0.2 andη¼0.1; the other known constants are given in Table 1.

E. Maréchal, O. Bastien / Journal of Theoretical Biology 361 (2014) 1–136

3.1.2. Asymptotic stability of steady states in the developed modelAs stated above, MCA being a first order sensitivity analysis in

the vicinity of a structurally stable fixed point, this fixed pointmust at least be linearly stable (Jones and Sleeman, 2003; Murray,2003). A steady state (us,vs) of the system (5) is said to be linearlystable if after a small perturbation, it returns to its original value.In mathematical terms, (us,vs) should respect this assumption ifthe linearized Jacobian J (or stability matrix J) for the steady statesatisfies tr(J)o0 and det(J)40, where tr(J) and det(J) are respec-tively the trace and the determinant of the matrix J (Jones andSleeman, 2003; Murray, 2003). In our model, we checked that allpossible steady states (us,vs) with positive values for us and vs wereindeed linearly stable (see Appendix A).

3.2. Estimates of the net incoming fluxes of PA and DAG at thestationary points in the inner envelope membrane, based onexperimental data

3.2.1. The net incoming flux of PAFollowing Dubots et al. (2010), we assumed that the PA

concentration in the IEM at stationary point [PA]s was reasonablyclose to, or in the range of the PA level measured in thylakoids,which was reported to be [PA]¼0.08 mol% (i.e., us¼0.013) (Fritz etal., 2007). This value is far below the estimated KM,PAP at 6 mol%.The level of DAG at the stationary point is usually considered to befar lower than 1 mol%, that is to say 0ovso1/6. A reasonablehypothesis for the DAG concentration in the IEM at stationarypoint is in the range of [DAG]¼0.1 mol%. Again this value is belowthe estimated KM,MGD1 at 1–5 mol%. The amount of enzymes istherefore unlikely to be limiting. Using the first equation in system(6), we can deduce the following inequalities: ðus=ð1þðb=6ÞÞþusÞo f oðus=1þusÞ.

With the following scenarios, we can therefore assess that,when [PA]s¼0.08 mol% and [DAG]s¼0.1 mol%, 0.0115o fo0.0132;when [PA]s¼0.08 mol% and [DAG]s¼1 mol%, 0.0055o fo0.0132;when [PA]s¼1 mol% and [DAG]s¼0.1 mol%, 0.1273o fo0.1429 andeventually when [PA]s¼1 mol% and [DAG]s¼1 mol%, 0.0642o fo0.1429. In dimensional terms, we have α¼ VM;PAPWf . Even ifW is unknown, it allows us to give an estimation of the incomingnet flux of PA because (i)W helps fixing a unit inconsistency to thisestimation (W is the number of MGD1 molecules divided by thenumber of lipids in the IEM) and (ii) W is not important for thesteady state value (W is related to the MGD1 enzyme quantity).So we have:

– when [PA]s¼0.08 mol% and [DAG]s¼0.1 mol%, 0.046oαo0.053;

– when [PA]s¼0.08 mol% and [DAG]s¼1 mol%, 0.022oαo0.053;– when [PA]s¼1 mol% and [DAG]s¼0.1 mol%, 0.509oαo0.572;– when [PA]s¼1 mol% and [DAG]s¼1 mol%, 0.257oαo0.572;

with α expressed in molecules of PA s�1 molecule of MGD1�1.In other words, if PA is present in the IEM at a concentration of

0.08 mol% or above, the PA influx can neither be null nor exceedingan upper limit. In the following representations, the PA flux wassimulated within this interval of lower and upper bounds. Atsteady state, the net influx of PA is also trivially equal to the flux of

DAG directly synthesized by the PAP. The value of the influx of PA,in molecules per time unit, is lower than the VM,MGD1, which isestimated about 20 molecules of DAG s�1.molecule of MGD1�1

(Table 1). This implies that a steady state can only exist if the influxis very small compared to the potential efficacy of the enzymes(estimated by the VM,MGD1, in Mendelian conditions) and that inthis steady state, the PAP/MGD1 enzymes efficiently utilize incom-ing PA without being limiting.

3.2.2. The net incoming flux of DAGUsing the second equation in the system (6), we can deduce

the following inequalities do ððg=6Þ=rþð1þð1=cusÞÞð1=6ÞÞ� fmin,where fmin is the lower bound of f deduced in the previous section.With KM,MGD1¼1 mol%, when [PA]s¼0.08 mol% and [DAG]s¼0.1 mol%, do0.1565; when [PA]s¼0.08 mol% and [DAG]s¼1 mol%, do0.2355; when [PA]s¼1 mol% and [DAG]s¼0.1 mol%,do0.2727 and eventually when [PA]s¼1 mol% and [DAG]s¼1 mol%, do1.3644.

In dimensional term, we have η¼VM,PAPWd, where W is stillunknown, giving for the DAG influx from the OEM:

– when [PA]s¼0.08 mol% and [DAG] s¼0.1 mol%, ηo0.58;– when [PA]s¼0.08 mol% and [DAG] s¼1 mol%, ηo0.84;– when [PA]s¼1 mol% and [DAG] s¼0.1 mol%, ηo0.97;– when [PA]s¼1 mol% and [DAG] s¼1 mol%, ηo4.87;

with η in molecules of DAG.s�1.molecules of MGD1�1.At steady state, these upper bound values of η are much higher

than the flux of DAG produced by the PAP, (equal to α, see above),suggesting that for an estimate of a DAG concentration in the IEM,most of it could be imported rather than synthesized in situ. The[DAG] value at 1 mol% is known to be an overestimate. With thisvalue, regardless the concentration of PA, the total influx of DAG isnevertheless below the turnover number of MGD1, i.e. 20 mole-cules of DAG s�1 molecule of MGD1�1 (Table 1), implying that asteady state can only exist if the influx is small compared to thepotential efficacy of the enzymes and that in this steady state, theMGD1 enzymes can also efficiently utilize incoming DAG, withoutbeing limiting.

3.3. Variation of the steady state in function of the incoming netfluxes of PA and DAG

Using the intervals for f and d obtained from Section 3.2, wecould compute all possible steady states for [PA] and [DAG] forvarious values of f and d in these intervals (Fig. 5). Steady statesvalues for PA (respectively DAG) define a 3-dimensional surface asa function of PA and DAG influxes in the system, i.e. α and η (Fig. 5,shown in panels A1, A2, A3, D1 and D2, respectively B1, B2, B3, C1ad C2). The shapes of [PA]S, and [DAG]S highlight clear differences.

Concerning [PA]S, when the incoming net flux of PA increases,the steady state concentration of PA increases nearly linearly.Likewise, when the incoming net flux of DAG increases, [PA]Svalue increases nearly linearly. The combination of DAG and PAinfluxes in the system has a simple additive effect. We can deducethe following features:

– [PA]S, is apparently not an important point of regulation of thesystem. It roughly follows the influx of PA as well as that ofDAG, as a consequence of the inhibition of PAP by this product.

– Although the system has the potential efficacy to convert PAinto DAG and then into MGDG, the presence of the tworegulatory loops prevents the IEM to be fully deprived in PA,when either PA or DAG incoming fluxes increase. This feature ofthe system implies an accumulation of PA to levels that are not

Table 2Necessary conditions for the existence of steady states in the developed model.All parameters are defined in the text.

Parameters Conditions Estimations (units)

α αoW.VM,PAP αo4 (PA mole fraction) s�1

αþη αþηoW.VM,MGD1 αþηo20 (PAþDAG mole fraction) s�1

E. Maréchal, O. Bastien / Journal of Theoretical Biology 361 (2014) 1–13 7

Fig. 5. Numerical simulations of the steady state values for [PA]S and [DAG]S in function of the incoming net fluxes of PA and DAG. Steady states values for [PA] (panels A1, A2,A3, D1 and D2) and [DAG] (panels B1, B2, B3, C1 ad C2) define a 3-dimensional surface as a function of PA and DAG influxes in the system, i.e. α and η in the model. The (α,η)values corresponding to a fixed value of [PA]S at 0.08 mol%, 0.5 mol% or 1.0 mol% define curves shown in black in (A1, A2, and A3 respectively) and the corresponding [DAG]Svalues define curves shown in white in (B1, B2 and B3 respectively). Vice versa, the (α,η) values corresponding to a fixed value of [DAG]S at 0.1 mol% or 1.0 mol% define curvesshown in white in (C1, and D2 respectively) and the corresponding [PA]S values define curves shown in black in (D1 and D2 respectively). Scale bars indicate mol% valueswith color codes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

E. Maréchal, O. Bastien / Journal of Theoretical Biology 361 (2014) 1–138

plausible. To reconcile numerical simulations with experimen-tal measures, one should hypothesize that PA should beexhausted from the IEM for another purpose than galactolipidsynthesis. Given the physiological conditions, the reasonablescenario of a [PA]S¼0.08 mol% in plants growing in phosphate-replete conditions (what is currently considered as a “normal”growth condition), defines an isocline curve (Fig. 5A1,B1),which is only consistent with very low influxes of PA(Fig. 5A1) and of DAG (Fig. 5B1). Given the proportion ofeukaryotic MGDG species, originating from an extraplastidialPA or DAG structure, the 0.08 mol% value for [PA]S could not beconsistent for a developing chloroplast. Simulations of [PA]S at0.5 mol% or even 1 mol% are compatible with increased influxesof PA and/or DAG (Fig. 5A2, A3), but in those cases, the influx ofPA should always be combined with an influx of DAG, acondition which is strictly necessary at [PA]S¼1 mol% implyingthat the influx of DAG can simply not be null (Fig. 5, B3).

Concerning [DAG]S, the shape of the surface is completelydifferent. If the incoming net flux of DAG is null, there is nearlyno possible accumulation of DAG, whatever the influx of PA. Whenthe influx of DAG increases, then [DAG]S increases. By contrast,when the influx of PA increases, although PA is the substrate forthe synthesis of DAG by the action of the PAP, the [DAG]S decreases.We can then distinguish two situations: either the influx of PA issmall and then [DAG]S can reach the 10 mol% range, or the influx ofPA is higher and then, DAG does not accumulate. Both DAG and PAinfluxes appear as necessary and they trigger therefore a dualeffect, based on the status of the PA influx. We can deduce thefollowing features:

– [DAG]S is clearly the parameter that is controlled in the system.It is mainly maintained to the lowest level as possible, exceptwhen PA influx is nearly shut down. This is consistent with PAbeing necessary to switch MGD1 on (Dubots et al., 2010). DAGinflux can then increase without triggering a strong accumula-tion of DAG in the IEM, as long as a PA influx is maintained inthe system. This feature is well adapted to a massive import ofDAG, for instance of eukaryotic origin, thus blocking the PAP,leaving PA for other purposes, and directing DAG toward MGDGsynthesis.

– A simulation of [DAG]S values at 1 mol% is only consistent witha combined entry of both PA and DAG in the system (Fig. 5, C2)and it implies an accumulation of PAwithin the IEM (Fig. 5, D2).This accumulation observed in the model can only be recon-ciled with the analysis of PA levels in chloroplasts in physiolo-gical conditions as long as an alternative utilization of PA ishypothesized, this alternative pathway being simply thebiosynthesis of PG.

3.4. Control coefficients

3.4.1. Flux control coefficientsThe computing of the flux control coefficients can be performed

analytically because kinetic expressions for MGD1 and PAP areknown from experiments and hence elasticity coefficients can becalculated (see Appendix A). Let Φ1, Φ2 and Φ3 variables corre-sponding to the fluxes via (1) α and η and β. Taking into accountthat the production of PA by the prokaryotic pathway wascombined with an import of PA from the eukaryotic pathway,which is fueled by driving forces outside the plastid, i.e. PLDs andPLCs, we assumed here that PA and DAG incoming net fluxes wereindependent of the PA and DAG IEM concentrations (i.e. elasticitycoefficients for α and η are null). Therefore, using the connectivitytheorem, the summation theorem together with the branching

point theorem (see Apendix), it can be demonstrated thatCΦ3PAP ¼ CΦ3

MGD1 ¼ 0, that is to say that the MGDG fluxes in the IMEdo not depend on the activity of the enzyme PAP and MGD1, andthat CΦ3

α ¼ a andCΦ3η ¼ 1�a where a is the fraction of the MGDG

fluxes in IME coming from α.

3.4.2. Concentration control coefficientsEven if the fluxes do not depend on the activity of PAP and

MGD1, the steady states for PA and DAG concentration coulddepend on the activity of these two enzymes. Using summationstheorems, connectivity theorems, branching point theorems andthe assumption that both PA and DAG incoming net fluxes areindependent of the PA and DAG inner membrane concentrations,we can calculate the eight concentration control coefficients forvarious PA and DAG incoming net fluxes (see Appendix A)

(i) Concentration control coefficients for α and ηː

CPAα ¼ εMGD1

DAG �aεPAPDAG

Ω;CPA

η ¼ εPAPDAGða�1ÞΩ

;

CDAGα ¼ aεPAPPA �εMGD1

PA

Ωand CDAG

η ¼ �εPAPPA ð1�aÞΩ

; ð8aÞ

(ii) Concentration control coefficients for PAP and DAGː

CPAPAP ¼ �εMGD1

DAG

Ω;CPA

MGD1 ¼εPAPDAG

Ω;CDAG

PAP ¼ εMGD1PA

Ω;

and CDAGMGD1 ¼ �εPAPPA

Ω; ð8bÞ

with Ω¼ εPAPPA εMGD1DAG �εPAPDAGε

MGD1PA and a is the fraction of the

MGDG outflow coming from α (Fig. 6).

Consistently with the variation of the steady state in function ofthe incoming net fluxes of PA and DAG described above, we thusobserved that both MGD1 and PAP activity could have an influenceon the steady state concentration of PA and DAG. Concerning thecontrol of PA by the enzymes of the system, both PAP and MGD1exert an effect on this precursor (Fig. 6, A and B), with negativecoefficient values, whatever the PA or DAG influxes, and with asub-planar shape indicative of the existence of a simple mode ofcontrol. PAP and MGD1 also exert a control on DAG (Fig. 6, Cand D), with positive and negative coefficient values respectively.The shapes of the corresponding surfaces are strikingly different,and show multiple modes of control, especially at low values of PA(peak in Fig. 6, C and well in Fig. 6, D). These reflect the role of PA asa switch in the system, acting more specifically on DAG.

It should be noticed that if activation of MGD1 by PA did notexist, we would have εMGD1

PA ¼ 0, and hence CDAGPAP ¼ 0. The activation

of MGD1 by PA is therefore critical in this system as it allows acontrol of the DAG concentration by PAP activity and a greaterindependence of the concentration control coefficients in respectto the magnitude of the α and η fluxes.

4. Conclusion

In conclusion, the model we constructed implies the followingdynamic properties in the biosynthetic pathway of MGDG in theIEM:

(1) Although we considered conditions in which, on the one hand,steady states concentrations of PA ([PA]S) and DAG ([DAG]S)were lower than the KM values of PAP and MGD1 respectively,and on the other hand, influxes of PA (α) and DAG (η) werelower than the kCAT of PAP and MGD1 (expressed per mole-cules of MGD1, see Table 1), the presence of the two regulatory

E. Maréchal, O. Bastien / Journal of Theoretical Biology 361 (2014) 1–13 9

loops (inhibition of PAP by DAG and stimulation of MGD1 byPA) allows an accumulation of both precursors.

(2) If the influx of PA increases, then the steady state concentra-tion of PA increases, mostly as a consequence of PAP retro-inhibition by its product. This increase of [PA]S, to levels thathave never been measured experimentally, occurs whateverthe origin of PA, should it be synthesized within the stroma viathe prokaryotic pathway or imported from the ER via theeukaryotic pathway. It is possible to reconcile the model withthe physiological context only by considering the existence ofan alternative pathway utilizing PA in the IEM, which istrivially the biosynthesis of PG. In Arabidopsis chloroplast PGis mostly of prokaryotic type. Nevertheless, it has been shownthat the PG biosynthesis pathway could tolerate eukaryoticsubstrates (Fritz et al., 2007). Therefore the low level ofeukaryotic-type PG in chloroplasts does not depend on thespecificity of PG-synthesizing enzymes. Thus, in silico numer-ical simulations presented here and experimental character-ization of PG molecular structure in vivo imply that the influxof PA should be mainly of prokaryotic type. A quantitativeimport of eukaryotic PA in the IEM is not supported bythe model.

(3) If the influx of DAG increases then [PA]S also increases, also dueto the retroinhibition of PAP by DAG. This is consistent with anorientation of PA toward an alternative utilization, i.e. PG,whereas DAG is guided toward galactolipid synthesis.

(4) The complete model seems to be simply designed to preventDAG accumulation, regardless of the magnitudes of the influ-xes of PA and DAG. This remarkable property is striking inFig. 5, in which [DAG]S is close to a flat plane except as soon as

α reaches a minimum threshold value. A very low influx of PAis thus necessary to turn the system on, by stimulating MGD1.An accumulation of DAG when the galactolipid:galactolipidgalactosyltransferase (or SFR2) is active for instance (Dorne etal., 1982; – now known to correspond to SFR2, Moellering etal., 2010) would then be possible by a complementary arrest ofPA influx. Vice versa, following SFR2 action, the IEM composi-tion could return to its initial situation by stimulating an influxof PA that could switch the PAP/MGD1 system on.

Taken together, the model we built based on experimental datais consistent with a preservation of PA produced by the prokar-yotic pathway for PG synthesis, an unlikely massive influx ofeukaryotic PA, and an influx of eukaryotic DAG in the IEM guidedtoward MGDG synthesis. This is also consistent with the loss of theprokaryotic MGDG synthesis reported in the evolution of Angios-perms after independent mutations that have occurred in numer-ous clades (Mongrand et al., 1998; Petroutsos et al., 2014). The TGDcomplex, made of the TGD1, TGD2 and TGD3 IEM subunits and aTGD4 component associated to the OEM has been proposed toimport eukaryotic PA in the chloroplast for galactolipid synthesis(Xu et al., 2010; Wang et al., 2012; Wang et al., 2013). Homologs ofall these subunits exist in the genome of plants that do notsynthesize prokaryotic galactolipids, i.e. lacking a conversion ofPA into DAG in the IEM (see Supplementary Fig. 1). Based onexperimental data and reasonable assumptions, the system wedescribe is indeed ready to “lose” the PAP activity upstreamgalactolipids, without impairing the synthesis of PG, using prokar-yotic PA, and the synthesis of MGDG, using eukaryotic DAG.

Fig. 6. Concentration control coefficients. In the MCA framework and using summations theorems, connectivity theorems, branching point theorems and the assumptionthat both PA and DAG incoming net fluxes are independent of the PA and DAG concentrations in the IEM, concentration control coefficients were computed for various PAand DAG incoming net fluxes Both MGD1 and PAP activity have an influence on the steady state concentration of PA and DAG. (A) Control of [PA]S by the activity of PAP.(B) Control of [PA]S by the activity of MGD1. (C) Control of [DAG]S by the activity of PAP. (D) Control of [DAG]S by the activity of MGD1.

E. Maréchal, O. Bastien / Journal of Theoretical Biology 361 (2014) 1–1310

Our prediction is fully consistent with in vitro analyses of PAturnover performed three decades ago on chloroplasts isolatedfrom 16:3- and 18:3 plants, labeled with [14C]-acetate (Gardinerand Roughan, 1983). In 16:3-chloroplasts, following fatty acidlabeling, PA is synthesized and rapidly hydrolyzed. By contrast,PA is rather stable in 18:3-chloroplasts, Chloroplasts isolated from18:3-plants do not form labeled DAG when incubated underappropriate conditions, whereas chloroplasts from 16:3-plantsdo (Gardiner et al., 1984). The molecular characterization of thePAP enzyme involved in this pathway would help investigating the16:3-to-18:3-plant transitions. Three plastid proteins, called LPPγ,LPPε1 and LPPε2, have been identified in Arabidopsis based ontheir homology with cyanobacterial lipid phosphatase-likesequences, and proposed to act as PAP enzymes (Nakamura etal., 2007), since they could partly complement yeast mutantsimpaired in PAP synthesis. Inconsistently with their possible rolein thylakoid lipid synthesis, loss-of-function of LPPε1/LPPε2 oroverexpression of all three genes did not lead to any visiblechange in lipid profiles (Nakamura et al., 2007). Homozygousknock out of LPPγ could not be obtained by Nakamura et al.(2007), apparently because of a failure of pollen tube growth.LPPγ might therefore play other function, which should make thissequence indispensable. Based on this lethal phenotype, it wastherefore proposed that LPPγ might be the IEM PAP. However, LPPγhomologs are clearly present in 18:3 plants such as Glycine max(XP_003537269) or Zea mais (NP_001146731), with a very highlevel of sequence similarity, which makes LPPγ an unlikelycandidate. Multiple homozygous Arabidopsis lines with a T-DNAinsertion in an LPPγ exon (SALK_055510C, SALK_048788C andSALK_005888C) seem to have been produced independently sincethe work by Nakamura et al. (2007), so the LPPγ precise functionmight be reexamined experimentally. More work is thereforeneeded to unambiguously define the enzyme that converts PAinto DAG in Arabidopsis chloroplasts.

Eventually, the role of TGD, if involved in eukaryotic PA importwould then be on the one hand to provide this lipid as a signalingmolecule to turn MGD1 on, and on the other hand to import otherlipid components. This conclusion revitalizes a long standingquestion, which is the very likely, but still undemonstrated, influxof DAG of eukaryotic origin in the IEM, either by direct import ofDAG in the IEM or conversion of eukaryotic PA or phospholipids inthe OEM prior transfer of eukaryotic DAG to the IEM.

Author disclosure statement

No competing financial interests exist.

Acknowledgments

We thank M.A. Block, L. Boudière, N. Glade and J. Jouhet forhelpful discussions. The authors were supported by the AgenceNationale de la Recherche (ANR-10-BLAN-1524, ReGal; ANR-12-BIME-0005, DiaDomOil), the Labex GRAL (Grenoble Alliance forIntegrated Structural Cell Biology) and Institut National de laRecherche Agronomique.

Appendix A

A.1. Analytic solution for the existence of a stationary point followingthe system of Eq. (6)

(i) The first Eq. in (6) corresponds to the resolution of Γ(u,v)¼0by solving bv¼ ðð1=f Þ�1Þu�1¼ Xu�1; where X ¼ ð1=f Þ�1:

(ii) The second of Eq. in (6) can now be reformulated taking intoaccount (i):

Δðu; vÞ ¼ 0 ) uð1þbvÞþu

� gvrþð1þð1=cuÞÞv þd¼ 0

) brXu�1

þ 1cu

þ1� gdþð1=Xþ1Þ ¼ 0

) brXu�1

þ 1cu

þY ¼ 0; where Y ¼ 1� gdþð1=Xþ1Þ

This leads to the following quadratic equation: XYcu2þ(brcþX–Yc)u–1¼0. Remarkably, the discriminant is always positive ifconditions for a stationary point are filled (see main text, statingthat X40 and Yo0). Then, the relevant solution for u and v aregiven by:

us ¼�ðbrcþX�YcÞ�

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiðbrcþX�YcÞ2þ4XYc

q2XYc

and vs ¼1b

1f�1

� �us�1

� �:

A.2. Linear stability of steady states in the developed model

The steady state (us,vs) is linearly stable if the linearizedJacobian J (or stability matrix J) for the steady state satisfies tr(J)o0 and det(J)40, where tr(J) and det(J) are respectively the traceand the determinant of the matrix J. With Γ(u,v) and Δ(u,v)defined in the model described in the main text, the Jacobianmatrix is:

J ¼∂Γ∂u

ðus; vsÞ ∂Γ∂v

ðus; vsÞ∂Δ∂u

ðus; vsÞ∂Δ∂v

ðus; vsÞ

0BB@

1CCA

where

∂Γ∂u

ðus; vsÞ ¼ Γu ¼ � 1þbvsð1þbvsþusÞ2

¼ �Γ1

∂Γ∂v

ðus; vsÞ ¼ Γv ¼ bus

ð1þbvsþusÞ2¼ Γ2

∂Λ∂u

ðus; vsÞ ¼ Δu ¼1þbvs

ð1þbvsþusÞ2� gv2scu2

s ðrþð1þð1=cusÞÞvsÞ2¼ Γ1�Δ1

∂Λ∂v

ðus; vsÞ ¼ Δv ¼ � bus

ð1þbvsþusÞ2� gr

ðrþð1þð1=cusÞÞvsÞ2¼ �Γ2�Δ2

The computing of tr(J) and det(J) with the previous notationgive:

trðJÞ ¼ ΓuþΔvo0

detðJÞ ¼ ΓuΔv�ΓvΔu ¼ �Γ1ð�Γ2�Δ2Þ�Γ2ðΓ1�Δ1Þ ¼ Γ1Δ2þΓ2Δ140

A.3. Elasticity coefficients

Elasticity Coefficients ευiAlare defined as the relative changes of a

reaction rate υi of the i-th reaction in response to a relative changeof a network species Al, that is to say ευiAl

¼ ð∂ log υi=∂ log AlÞ. In the

model illustrated in Fig. 3, both α and η are assumed to beindependent of the concentration of the species [PA] and [DAG],as a consequence, we haveεαPA ¼ εαDAG ¼ εηPA ¼ εηDAG ¼ 0. For other

elasticity coefficients, we have εPAPPA ¼ ð1þbv=ð1þbvþuÞÞ, εPAPDAG ¼�ðbv=ð1þbvþuÞÞ, εMGD1

PA ¼ vcurþð1þð1=cuÞÞv and εMGD1

DAG ¼ rrþð1þð1=cuÞÞv.

E. Maréchal, O. Bastien / Journal of Theoretical Biology 361 (2014) 1–13 11

A.4. Determination of the flux and concentration control coefficients

Let Φ1, Φ 2 and Φ3 be the variables corresponding to the steadystate fluxes via α, η and β. By the law of conservation of mass, atsteady state, the fluxes in each limb will be governed by therelationship: Φ1þΦ2 ¼Φ3. For the branched system we canexpress the terms of the summation and connectivity theoremswith respect to each flux. Here, we are mainly interested in theproduction of MGDG, i.e. flux Φ3 ¼ β. Hence, we have thefollowing three equations with four unknown flux controlcoefficients:

CΦ3α þCΦ3

PAPþCΦ3η þCΦ3

MGD1 ¼ 1

CΦ3α ευiPAþCΦ3

PAPευiPAþCΦ3

η ευiPAþCΦ3MGD1ε

υiPA ¼ 0

CΦ3α ευiDAGþCΦ3

PAPευiDAGþCΦ3

η ευiDAGþCΦ3MGD1ε

υiDAG ¼ 0

To solve the expression of CΦυi, we need another equation. The

same happens with the two sets of concentration controlcoefficients, where the summation theorem states that

∑iCAnυi

¼ 0 for both PA and DAG and the connectivity theorem

for concentrations states that ∑iCAnυiευiAm

¼ �1, if n¼m or

∑iCAnυiευiAm

¼ 0 otherwise. Again, for each metabolite, we havethree equations but four unknown flux control coefficients. To

solve the expressions of CAnυi, additional equations are needed.

Following the approach of Kacser (1983), Fell and Sauro (1985)and Sauro et al. (1987), we can derive a branching point theoremproviding the missing relations.

Starting from the pathway depicted in Fig. 3, let the fraction of Φ3

originating fromΦ1 be given by a¼Φ1/Φ3 and that originating fromΦ2

be given by 1–a¼Φ2/Φ3. We can consider a perturbation of thesystem, in which the concentration of the enzyme and/or transporterinvolved in the α process and PAP increases by a given value, implyingan increase of Φ1 and [DAG]. As a consequence, this perturbationcauses (in general) a decrease inΦ2 (relief of product inhibition, if any)and an increase in Φ3. We can restore the change in Φ3 by decreasingΦ2 such that [DAG] is restored to its pre-perturbation state. At the endof this perturbation/restoration process, Δ[DAG]¼0. Since we have notchanged MGD1 concentration, ΔΦ3¼0. In this scenario, the systemequations can be formulated as follows:

dΦ3 ¼∂Φ3

∂αΔαþ ∂Φ3

∂PAPΔPAPþ∂Φ3

∂βΔβ

Using the relation CΦ3υi

¼ ð∂Φ3=∂υiÞðυi=Φ3Þ, we have ðdΦ3=Φ3Þ ¼CΦ3α ðΔα=αÞþCΦ3

PAPðΔPAP=PAPÞþCΦ3η ðΔη=ηÞ.

Local equations can be expressed as ðΔΦ1=Φ1Þ ¼ ðΔα=αÞ ¼ðΔPAP=PAPÞ and ðdΦ2=Φ2Þ ¼ ðΔη=ηÞ, with ΔΦ1þΔΦ2 ¼ 0. Combiningthe system and the local equations and remembering that ΔΦ3 ¼ 0,we have:

dΦ3

Φ3¼ ðCΦ3

α þCΦ3PAPÞð1�aÞ�CΦ3

η a¼ 0:

This formula is a branching point theorem for our pathway.This expression is highly similar with the generalized formulagiven by Sauro et al. (1987). It is straightforward to develop thesame reasoning to provide expressions of branching point theo-rems for both DAG and PA concentration control coefficients.

Summation theorems, connectivity theorems and branchingpoint theorems for Φ3, PA and DAG can be encompassed

in a matrix formalism as demonstrated by Sauro et al. (1987).The equation to be solved is FG¼H where

F ¼

1 1 1 1εαPA εPAPPA εηPA εMGD1

PA

εαDAG εPAPDAG εηDAG εMGD1DAG

1�a 1�a �a 0

0BBB@

1CCCA, G¼

CΦ3α CPA

α CDAGα

CΦ3PAP CPA

PAP CDAGPAP

CΦ3η CPA

η CDAGη

CΦ3MGD1 CPA

MGD1 CDAGMGD1

0BBBBB@

1CCCCCA

and G¼

1 0 00 �1 00 0 �10 0 0

0BBB@

1CCCA.

This linear system can be solved directly by Cramer's rule(Habgood and Arel, 2012) to give the results given in eq. (8a,b).

Appendix B. Supporting information

Supplementary data associated with this article can be found inthe online version at http://dx.doi.org/10.1016/j.jtbi.2014.07.013.

References

Andersson, M.X., Kjellberg, J.M., Sandelius, A.S., 2004. The involvement of cytosoliclipases in converting phosphatidyl choline to substrate for galactolipid synth-esis in the chloroplast envelope. Biochim. Biophys. Acta 1684 (1-3), 46–53.

Awai, K., Maréchal, E., Block, M.A., Brun, D., Masuda, T., Shimada, H., Takamiya, K.,Ohta, H., Joyard, J., 2001. Two types of MGDG synthase genes, found widely inboth 16:3 and 18:3 plants, differentially mediate galactolipid syntheses inphotosynthetic and nonphotosynthetic tissues in Arabidopsis thaliana. Proc.Natl. Acad. Sci. USA 98 (19), 10960–10965.

Benning, C., Ohta, H., 2005. Three enzyme systems for galactoglycerolipid biosynth-esis are coordinately regulated in plants. J. Biol. Chem. 280 (4), 2397–2400.

Boudière, L., Botté, C.Y., Saidani, N., Lajoie, M., Marion, J., Bréhélin, L., Yamaryo-Botté, Y., Satiat-Jeunemaître, B., Breton, C., Girard-Egrot, A., Bastien, O., Jouhet,J., Falconet, D., Block, M.A., Maréchal, E., 2012. Galvestine-1, a novel chemicalprobe for the study of the glycerolipid homeostasis system in plant cells. Mol.Biosyst. 8 (8), 2023–2035.

Boudière, L., Michaud, M., Petroutsos, D., Rébeillé, F., Falconet, D., Bastien, O., Roy, S.,Finazzi, G., Rolland, N., Jouhet, J., Block, M.A., Maréchal, E., 2014. Glycerolipids inphotosynthesis: Composition, synthesis and trafficking. Biochim. Biophys. Acta1837 (4), 470–480.

Block, M.A., Douce, R., Joyard, J., Rolland, N., 2007. Chloroplast envelope mem-branes: a dynamic interface between plastids and the cytosol. Photosynth. Res.92 (2), 225–244.

Cohen, E.R., 1996. The Physics Quick Reference Guide. AIP Press, New-York.Covés, J., Joyard, J., Douce, R., 1988. Lipid requirement and kinetic studies of

solubilized UDP-galactose:diacylglycerol galactosyltransferase activity fromspinach chloroplast envelope membranes. Proc. Natl. Acad. Sci. USA 85 (14),4966–4970.

Cornish-Bowden, A., Cárdenas, M.L., 1990. Control of Metabolic Processes. PlenumPress, New-York.

Cornish-Bowden, A., 2004. Fundamentals of Enzyme Kinetics, third ed. PortlandPress, London.

Cruz-Ramírez, A., Oropeza-Aburto, A., Razo-Hernández, F., Ramírez-Chávez, E.,Herrera-Estrella, L., 2006. Phospholipase DZ2 plays an important role inextraplastidic galactolipid biosynthesis and phosphate recycling in Arabidopsisroots. Proc. Natl. Acad. Sci. USA 103 (17), 6765–6770.

Dantzig, G.B., 1998. Linear programming and extensions. Princeton University Press,Princeton.

Dennis, E.A., Darke, P.L., Deems, R.A., Kensil, C.R., Plückthun, A., 1981. Cobra venomphospholipase A2: a review of its action toward lipid/water interfaces. Mol. CellBiochem. 36 (1), 37–45.

Dorne, A.J., Block, M.A., Joyard, J., Douce, R., 1982. The galactolipid: galactolipidgalactosyltransferase is located on the outer surface of the outer membrane ofthe chloroplast envelope. FEBS Lett. 145 (1), 30–34.

Dubots, E., Audry, M., Yamaryo, Y., Bastien, O., Ohta, H., Breton, C., Maréchal, E.,Block, M.A., 2010. Activation of the chloroplast monogalactosyldiacylglycerolsynthase mgd1 by phosphatidic acid and phosphatidylglycerol. J. Biol. Chem.285 (9), 6003–6011.

Fell, D.A., Sauro, H.M., 1985. Metabolic control and its analysis. Eur. J. Biochem. 148(3), 555–561.

Fell, D.A., 2005. Enzymes, metabolites and fluxes. J. Exp. Bot. 56 (410), 267–272.Frentzen, M., Heinz, E., McKeon, T.A., Stumpf, P.K., 1983. Specificities and selectiv-

ities of glycerol-3-phosphate acyltransferase and monoacylglycerol-3-phosphate acyltransferase from pea and spinach chloroplasts. Eur. J. Biochem.129 (3), 629–636.

E. Maréchal, O. Bastien / Journal of Theoretical Biology 361 (2014) 1–1312

Fritz, M., Lokstein, H., Hackenberg, D., Welti, R., Roth, M., Zähringer, U., Fulda, M.,Hellmeyer, W., Ott, C., Wolter, F.P., Heinz, E., 2007. Channeling of eukaryoticdiacylglycerol into the biosynthesis of plastidial phosphatidylglycerol. J. Biol.Chem. 282 (7), 4613–4625.

Gardiner, S.E., Heinz, E., Roughan, P.G., 1984. Rates and products of long-chain fattyacid synthesis from [14-C]acetate in chloroplasts isolated from leaves of 16:3and 18:3 plants. Plant Physiol. 74 (4), 890–896.

Gardiner, S.E., Roughan, P.G., 1983. Relationship between fatty-acyl composition ofdiacylgalactosylglycerol and turnover of chloroplast phosphatidate. Biochem. J.210 (3), 949–952.

Habgood, K., Arel, I., 2012. A condensation-based application of Cramer's rule forsolving large-scale linear systems. J. Discrete Algorithms 10, 98–109.

Heinrich, R., Rapoport, T.A., 1974. A linear steady-state treatment of enzymaticchains. General properties, control and effector strength. Eur. JBiochem. 42 (1),89–95.

Heinrich, R., Schuster, S., 1996. The Regulation of Cellular Systems. Chapman andHall, New York.

Heinz, E., Roughan, P.G., 1983. Similarities and differences in lipid metabolism ofchloroplasts isolated from 18:3 and 16:3 plants. Plant Physiol. 72 (2), 273–279.

Higgins, J., 1963. Analysis of sequential reactions. Ann. N.Y. Acad. Sci. 108, 305–321.Jarvis, P., Dörmann, P., Peto, C.A., Lutes, J., Benning, C., Chory, J., 2000. Galactolipid

deficiency and abnormal chloroplast development in the arabidopsis mgdsynthase 1 mutant. Proc. Natl. Acad. Sci. USA 97 (14), 8175–8179.

Jones, D.S., Sleeman, B.D., 2003. Differential Equations and Mathematical Biology.Chapman & Hall/CRC, New-York.

Kacser, H., Burns, J.A., 1973. The control of flux. Symp. Soc. Exp. Biol. 27, 65–104.Kacser, H., 1983. The control of enzyme systems in vivo: elasticity analysis of the

steady state. Biochem. Soc. Trans. 11 (1), 35–40.Kobayashi, K., Kondo, M., Fukuda, H., Nishimura, M., Ohta, H., 2007. Galactolipid

synthesis in chloroplast inner envelope is essential for proper thylakoidbiogenesis, photosynthesis, and embryogenesis. Proc. Natl. Acad. Sci. USA 104(43), 17216–17221.

Malherbe, A., Block, M.A., Joyard, J., Douce, R., 1992. Feedback inhibition ofphosphatidate phosphatase from spinach chloroplast envelope membranes bydiacylglycerol. J. Biol. Chem. 267 (33), 23546–23553.

Maréchal, E., Block, M.A., Joyard, J., Douce, R., 1994a. Kinetic properties ofmonogalactosyldiacylglycerol synthase from spinach chloroplast envelopemembranes. J. Biol. Chem. 269 (8), 5788–5798.

Maréchal, E., Block, M.A., Joyard, J., Douce, R., 1994b. Comparison of the kineticproperties of MGDG synthase in mixed micelles and in envelope membranesfrom spinach chloroplast. FEBS Lett. 352 (3), 307–310.

Michaelis, L., Menten, M.L., 1913. Die Kinetik der Invertinwirkung. Biochem. Z 49,333–369.

Miège, C., Maréchal, E., Shimojima, M., Awai, K., Block, M.A., Ohta, H., Takamiya, K.,Douce, R., Joyard, J., 1999. Biochemical and topological properties of type AMGDG synthase, a spinach chloroplast envelope enzyme catalyzing the synth-esis of both prokaryotic and eukaryotic MGDG. Eur. J. Biochem. 265 (3),990–1001.

Mongrand, S., Bessoule, J.J., Cabantous, F., Cassagne, C., 1998. The C16: 3/C18:3 fattyacid balance in photosynthetic tissues from 468 plant species. Phytochemistry49 (4), 1049–1064.

Mongrand, S., Cassagne, C., Bessoule, J.J., 2000. Import of lyso-phosphatidylcholineinto chloroplasts likely at the origin of eukaryotic plastidial lipids. Plant Physiol.122 (3), 845–852.

Moellering, E.R., Muthan, B., Benning, C., 2010. Freezing tolerance in plants requireslipid remodeling at the outer chloroplast membrane. Science 330 (6001),226–228.

Murray, J.D., 2003. Mathematical Biology. I: An Introduction, Third ed. Springer,Berlin.

Nakamura, Y., Awai, K., Masuda, T., Yoshioka, Y., Takamiya, K., Ohta, H., 2005.A novel phosphatidylcholine-hydrolyzing phospholipase C induced by phos-phate starvation in Arabidopsis. J. Biol. Chem. 280 (9), 7469–7476.

Nakamura, Y., Tsuchiya, M., Ohta, H., 2007. Plastidic phosphatidic acid phosphatasesidentified in a distinct subfamily of lipid phosphate phosphatases withprokaryotic origin. J. Biol. Chem. 282 (39), 29013–29021.

Nakamura, Y., Koizumi, R., Shui, G., Shimojima, M., Wenk, M.R., Ito, T., Ohta, H.,2009. Arabidopsis lipins mediate eukaryotic pathway of lipid metabolism andcope critically with phosphate starvation. Proc. Natl. Acad. Sci. USA 106 (49),20978–20983.

Nishiyama, Y., Hardré-Liénard, H., Miras, S., Miège, C., Block, M.A., Revah, F., Joyard, J.,Maréchal, E., 2003. Refolding from denatured inclusion bodies, purification tohomogeneity and simplified assay of MGDG synthases from land plants. ProteinExpr. Purif. 31 (1), 79–87.

Oropeza-Aburto, A., Cruz-Ramírez, A., Acevedo-Hernández, G.J., Pérez-Torres, C.A.,Caballero-Pérez, J., Herrera-Estrella, L., 2012. Functional analysis of the Arabi-dopsis PLDZ2 promoter reveals an evolutionarily conserved low-Pi-responsivetranscriptional enhancer element. J. Exp. Bot. 63 (5), 2189–2202.

Peano, G., 1889. Sur le déterminant Wronskien. Mathesis 9, 75–76.Petroutsos, D., Amiar, S., Abida, H., Dolch, L.J., Bastien, O., Rébeillé, F., Jouhet, J.,

Falconet, D., Block, M.A., McFadden, G.I., Bowler, C., Botté, C., Maréchal, E., 2014.Evolution of galactoglycerolipid biosynthetic pathways – from cyanobacteria toprimary plastids and from primary to secondary plastids. Prog. Lipid Res. 54C,68–85.

Salinas, D.G., Reyes, J.G., De la Fuente, M., 2011. Lipid metabolizing enzymeactivities modulated by phospholipid substrate lateral distribution. Bull. Math.Biol. 73 (9), 2045–2067.

Sauro, H.M., Small, J.R., Fell, D.A., 1987. Metabolic control and its analysis. Eur. J.Biochem. 165 (1), 215–221.

Segel, L.A., 1972. Simplication and scaling. SIAM Rev. 14, 547–571.The MathWorks, 2013. MATLAB and Optimization Toolbox Release. Inc, 2013a.

Natick, Massachusetts, United States.Wang, Z., Xu, C., Benning, C., 2012. TGD4 involved in endoplasmic reticulum-to-

chloroplast lipid trafficking is a phosphatidic acid binding protein. Plant J. 70(4), 614–623.

Wang, Z., Anderson, N.S., Benning, C., 2013. The phosphatidic acid binding site ofthe Arabidopsis trigalactosyldiacylglycerol 4 (TGD4) protein required for lipidimport into chloroplasts. J. Biol. Chem. 288 (7), 4763–4771.

Xu, C., Yu, B., Cornish, A.J., Froehlich, J.E., Benning, C., 2006. Phosphatidylglycerolbiosynthesis in chloroplasts of Arabidopsis mutants deficient in acyl-ACPglycerol-3- phosphate acyltransferase. Plant J. 47 (2), 296–309.

Xu, C., Moellering, E.R., Muthan, B., Fan, J., Benning, C., 2010. Lipid transportmediated by Arabidopsis TGD proteins is unidirectional from the endoplasmicreticulum to the plastid. Plant Cell Physiol. 51 (6), 1019–1028.

Yu, B., Xu, C., Benning, C., 2002. Arabidopsis disrupted in SQD2 encoding sulfolipidsynthase is impaired in phosphate-limited growth. Proc. Natl. Acad. Sci. U S A99 (8), 5732–5737.

E. Maréchal, O. Bastien / Journal of Theoretical Biology 361 (2014) 1–13 13