Computational inference of the structure and regulation of the lignin pathway in Panicum virgatum.

BACKGROUND: Switchgrass is a prime target for biofuel production from inedible plant parts and has been the subject of numerous investigations in recent years. Yet, one of the main obstacles to effective biofuel production remains to be the major problem of recalcitrance. Recalcitrance emerges in part from the 3-D structure of lignin as a polymer in the secondary cell wall. Lignin limits accessibility of the sugars in the cellulose and hemicellulose polymers to enzymes and ultimately decreases ethanol yield. Monolignols, the building blocks of lignin polymers, are synthesized in the cytosol and translocated to the plant cell wall, where they undergo polymerization. The biosynthetic pathway leading to monolignols in switchgrass is not completely known, and difficulties associated with in vivo measurements of these intermediates pose a challenge for a true understanding of the functioning of the pathway.
RESULTS: In this study, a systems biological modeling approach is used to address this challenge and to elucidate the structure and regulation of the lignin pathway through a computational characterization of alternate candidate topologies. The analysis is based on experimental data characterizing stem and tiller tissue of four transgenic lines (knock-downs of genes coding for key enzymes in the pathway) as well as wild-type switchgrass plants. These data consist of the observed content and composition of monolignols. The possibility of a G-lignin specific metabolic channel associated with the production and degradation of coniferaldehyde is examined, and the results support previous findings from another plant species. The computational analysis suggests regulatory mechanisms of product inhibition and enzyme competition, which are well known in biochemistry, but so far had not been reported in switchgrass. By including these mechanisms, the pathway model is able to represent all observations.
CONCLUSIONS: The results show that the presence of the coniferaldehyde channel is necessary and that product inhibition and competition over cinnamoyl-CoA-reductase (CCR1) are essential for matching the model to observed increases in H-lignin levels in 4-coumarate:CoA-ligase (4CL) knockdowns. Moreover, competition for 4-coumarate:CoA-ligase (4CL) is essential for matching the model to observed increases in the pathway metabolites in caffeic acid O-methyltransferase (COMT) knockdowns. As far as possible, the model was validated with independent data.

Background

About 440 million years ago plants started to leave the oceans and inhabit land [1, 2]. The emergence of n class="Chemical">lignin durinpan>g this time was anpan> adaptationpan> to the new enpan>vironpan>menpan>t anpan>d, specifically, a responpan>se to gravity anpan>d to limitationpan>s inpan> accessinpan>g pan> class="Chemical">water. The new life also demanded plants to store water and develop systems of water transfer. The plant furthermore needed to grow in height in order to have enough access to sunlight and oxygen. Plants ultimately accomplished these multiple tasks through their xylem structures, of which lignin is a key constituent. Lignin is a phenolic polymer that is woven around and between cellulose and hemicellulose within the secondary cell wall; it provides strength and facilitates water transfer in plants. A consequence of these significant benefits for plants is that lignin is very difficult to decompose, because it is an irregular polymer that contains aromatic rings. This resistance against decomposition and digestion is known as recalcitrance. It is arguably the most important barrier to industrializing second-generation biofuels, and in particular the production of ethanol from inedible plant parts as sustainable and affordable biofuels, because recalcitrance necessitates additional treatment steps, such as hot acid or ammonia baths, to loosen the lignin structure [3-5]. These steps require time and expense and therefore reduce feasibility and cost effectiveness. Moreover, most of the pretreatments are not environmentally friendly [6, 7]. Outside the biofuel industry, recalcitrance affects forage digestibility, and progress toward reducing recalcitrance could have a significant impact on the cattle and sheep industry [8].

Numerous attempts have been made in recent times to manipulate the n class="Chemical">lignin conpan>tenpan>t anpan>d compositionpan> inpan> canpan>didate planpan>ts for biofuel productionpan>. Manpan>y of these studies relied onpan> the assumptionpan> that the pan> class="Chemical">lignin biosynthesis pathway was known. However, this is not necessarily the case, especially in understudied plant species, and the precise pathway structure is often unclear and requires dedicated research for such species. For instance, Selaginella moellendorffi and Medicago truncatula have basically similar lignin pathways, which however differ in some of their metabolic branch points as well as their enzyme properties [9-11]. Beyond the topological structure, it is not surprising that different species have evolved distinct regulatory control patterns. The immediate consequence of such discrepancies for the biofuel industry is that the direct extrapolation of knowledge, methods and treatments from one species to another is not necessarily valid. Moreover, it is well known that pathway systems are highly nonlinear and difficult to predict with intuition alone. A feasible strategy is therefore to employ computational approaches of systems biology and metabolic engineering.

The design of suitable models for this purpose is not trivial. First, it is generically unclear which mathematical representations are optimal for describing a natural system. Second, one cannot be sure that information or data from one species can be assumed to be valid in another species, even if the two are closely related. Similarly, it has been shown many times that data obtained in vitro are not necessarily applicable in vivo [10-14]. At the same time, species-specific experiments are time consuming and expensive. Mechanistic models based on enzyme kinetics seem to be an intriguing choice, but it has been shown that mechanistic models are not always good solutions, for instance, if parameter values and enzymatic rate laws are based on strong assumptions like bulk reactivity that are not necessarily satisfied in vivo [12]. An alternative that was recently proposed is the characterization of in vivo-like kinetics [13], which however is costly and time consuming and would still require extensive validation, which however is seldom truly achieved [12]. An additional challenge for the design of models is the scarcity and quality of test and validation data, which pose a significant obstacle to all analyses of relatively understudied species.

In this study we analyze the n class="Chemical">lignin biosynpan>thesis pathway inpan> pan> class="Species">switchgrass, Panicum virgatum, with computational means of systems biology. The analysis is based on a dataset from stem and tiller tissue that consists of the lignin content (H, G and S lignin) and the S/G lignin ratio in wild type and in four transgenic lines (4CL, CCR1, CAD and COMT knockdowns). To some degree, details of the in vitro kinetics of some of the pathway enzymes have also been determined by one of our labs. Our approach here is to develop computational models that characterize the structure and regulatory control patterns of lignin biosynthesis in P. virgatum at a systemic level. The goals of this modeling approach are, first, to explain the experimental results from wild type and transgenic lines and, second, to devise a rational basis for strategies to manipulate the pathway toward reduced recalcitrance.

Results

The results are described in a sequence that follows our sten class="Chemical">p-by-step model design anpan>d conpan>veys our rationpan>ale for utilizinpan>g the observationpan>s to remediate discrepanpan>cies with the data anpan>d for suggestinpan>g the inpan>vestigationpan> of new features to the model inpan> the next step of the anpan>alysis. We beginpan> by assessinpan>g the pathway structure inpan> pan> class="Species">switchgrass as it is alleged in the current literature. Next, we examine possible channeling of CCR/CAD, which has been reported for the lignin pathway in alfalfa [5, 14], but not in switchgrass. Even accounting for the possibility of channeling, the experimental data regarding H lignin cannot be captured at this point. Thus, we investigate the effects of product inhibition and competitive inhibition. In the next phase, 4CL inhibition is added as a potential explanation for the accumulation of 4CL substrates, along with a simultaneous decrease in coniferaldehyde in the COMT knockdown. Finally, principal component analysis is performed to investigate the distribution of parameters within the high-dimensional parameter space and to reduce the feasible subspace of parameter values. The results section ends with a validation of the model.

Reaction system of lignin biosynthesis in switchgrass

The traditionally accepted n class="Chemical">lignin biosynpan>thesis pathway branpan>ches at pan> class="Chemical">p-coumaroyl CoA to provide S and G-lignin precursors (Fig. 1). The hexagon in this figure shows the details of this branch point. It was also previously assumed, based on studies in the dicots A. thaliana and N. benthamiana, that p-coumaroyl CoA is converted to p-coumaroyl shikimate and p-coumaroyl quinic acid by HCT. Subsequently, both products, p-coumaroyl shikimate and p-coumaroyl quinic acid, were shown to be converted to caffeoyl shikimate and caffeoyl quinic acid, respectively [15]. The enzyme for these unidirectional reactions is C3′H. Downstream, HCT was proposed to operate in the reverse direction to convert caffeoyl shikimate and caffeoyl quinic acid into caffeoyl-CoA.

Fig. 1

Lignin biosynthesis pathway. Dashed arrows represent the traditionally accepted pathway of lignin biosynthesis, while the arrow from caffeoyl shikimate to caffeic acid captures a newly discovered enzymatic activity [39] now known to be present in switchgrass. Caffeoyl shikimate esterase turns caffeoyl shikimate into caffeic acid and circumvents the previously accepted route. 4CL has recently been shown to exhibit activity towards caffeic acid and ferulic acid in switchgrass by which a new network topology is introduced for switchgrass lignin biosynthesis. Note that tyrosine is shown here, but not included in the model

n class="Chemical">Lignin biosynpan>thesis pathway. Dashed arrows represenpan>t the traditionpan>ally accepted pathway of pan> class="Chemical">lignin biosynthesis, while the arrow from caffeoyl shikimate to caffeic acid captures a newly discovered enzymatic activity [39] now known to be present in switchgrass. Caffeoyl shikimate esterase turns caffeoyl shikimate into caffeic acid and circumvents the previously accepted route. 4CL has recently been shown to exhibit activity towards caffeic acid and ferulic acid in switchgrass by which a new network topology is introduced for switchgrass lignin biosynthesis. Note that tyrosine is shown here, but not included in the model

A recent study demonstrated that this pathway organization is unlikely to occur in n class="Species">switchgrass [16]. Based onpan> kinpan>etic measuremenpan>ts of PvHCT1a, PvHCT2a anpan>d PvHCT-Like1, it was shownpan> that pan> class="Chemical">caffeoyl shikimate is not converted to caffeoyl-CoA by the reverse HCT reaction, but is more likely converted into caffeic acid through caffeoyl shikimate esterase, and that this step is actually the main route of mass transfer into the pathway towards S and G monolignols. As indicated with dashed arrows in Fig. 1, HCT is not active in the formation of caffeoyl-CoA. This new information helps us reduce the steps in Fig. 1. It has furthermore been suggested that cinnamic acid is a precursor for salicylic acid; this process is represented by the thick grey arrow [5]. Similarly, a considerable portion of ferulic acid leaves the pathway [17]. Finally, the efflux out of p-coumaric acid acts to avoid accumulation of the metabolite in the 4CL knockdown strain (Fig. 1). These simplifications yield the pathway diagram in Fig. 2.

Fig. 2

Revised and simplified pathway in switchgrass. By eliminating HCT from the diagram in Fig. 1 and adding CSE, the pathway system becomes simpler. The right branch in the grey box in Fig. 1 is merged into an efflux and the left branch is simplified to a one-step process. It is hypothesized that a specific functional channel could facilitate the conversion of feruloyl-CoA into coniferyl alcohol. Such a channel could be the result of co-localization of the involved pathway enzymes

Revised and simplified pathway in n class="Species">switchgrass. By eliminpan>atinpan>g HCT from the diagram inpan> Fig. 1 anpan>d addinpan>g pan> class="Gene">CSE, the pathway system becomes simpler. The right branch in the grey box in Fig. 1 is merged into an efflux and the left branch is simplified to a one-step process. It is hypothesized that a specific functional channel could facilitate the conversion of feruloyl-CoA into coniferyl alcohol. Such a channel could be the result of co-localization of the involved pathway enzymes

At this point, it is not entirely clear whether the n class="Chemical">lignin pathway inpan> pan> class="Species">switchgrass contains caffeyl aldehyde. It appears that this is not the case, and the following analysis assumes that caffeyl aldehyde is indeed not produced. Nonetheless, since other species do generate this intermediate, the Additional file 1: Text S1 analyzes this case.

Large-scale simulation studies with this pathway structure lead to irreconcilable differences between the experimental data and the model results, which indicate that the model has genuine flaws. In particular, the dynamics of the different n class="Chemical">lignin species cannot be explainpan>ed for the various transgenpan>ics (data not shownpan>).

Channeling

Experimental and theoretical work in alfalfa has suggested that functional enzymatic channeling likely occurs at the n class="Chemical">coniferaldehyde node [5, 14]. Accordinpan>g to this suggestionpan>, the “G-chanpan>nel” facilitates the use of pan> class="Chemical">feruloyl-CoA for the production of coniferyl alcohol, which is the precursor of the G monolignol (Fig. 2). We investigate the same channeling hypothesis here as a possibility. Specifically, we use pertinent experimental data from switchgrass to analyze the feasibility of different hypothetical pathway topologies. The potential existence of a functional complex consisting of CCR1/CAD leads to three possible pathway topologies that satisfy the requirement of mass conservation (Fig. 3).

Fig. 3

Topological Configurations. Three pathway structures are plausible when a CCR1/CAD channel is considered. Configuration 2 lacks the channel, while the other two configurations represent alternatives involving the channel

Each of these so-far unregulated topologies was modeled as a generalized mass action (GMA) model, whose parameter values were obtained with a sophisticated large-scale sampling scheme (see “Methods”). Although all topologies were found to be consistent with most of the experimental results, no topology was compatible with the accumulation of n class="Chemical">H lignin inpan> 4CL knpan>ockdownpan> tranpan>sgenpan>ics (Table 1); this situationpan> could not be simulated by anpan>y of the canpan>didate models, regardless of the presenpan>ce or absenpan>ce of the chanpan>nel. This stronpan>g result suggests the existenpan>ce of regulatory mechanpan>isms, anpan>d conpan>siderinpan>g the structure of the pathway anpan>d the branpan>ch toward pan> class="Chemical">H lignin in particular, we decided to analyze the possible role of product inhibition, which is frequently found in pathway systems in vivo.

Table 1

Fold change in lignin monomers, total lignin, and S/G in transgenic plants relative to wild-type plants

	4CL knockdown 40 % [3]	CCR knockdown 50 % [40]	COMT knockdown 30 % [5]	CAD knockdown 30 % [41]
Down-regulation	27–95 %	Up to 75 %	Up to 90 %	55–86 %
H lignin	1.82	NR	NR	NR
G lignin	0.53	~0.75	0.76–0.98	0.67–0.83
S lignin	1.00	~0.75	0.42–0.96	0.58–0.87
Total lignin	0.78	~0.75	0.84–0.96	0.78–0.86
S/G	Increased	Increased	Decreased	Decreased

NR not reported

Product inhibition

Experimental results from transgenic plants have demonstrated that n class="Chemical">H lignin accumulates whenpan> the enpan>zyme 4CL is downpan>-regulated [3]. Anpan>alyzinpan>g this inpan>itially counpan>terinpan>tuitive observationpan> closer suggests that there might be a wave of accumulationpan> inpan> the metabolites precedinpan>g pan> class="Chemical">H lignin. Such a wave can be explained with product inhibition (Fig. 4). When an enzyme is down-regulated, the corresponding substrate accumulates. The secondary effect is that the accumulated substrate is by itself a product of a previous reaction whose increased concentration decreases its own rate of production. This backward cascade has an upstream domino effect along the pathway and, depending on the kinetics of the reactions, can lead to the accumulation of upstream metabolites. This observation can be explained by the following chain of events: Down-regulating 4CL leads to a decrease in the products of this enzyme, i.e., p-coumaroyl-CoA, caffeoyl-CoA, and feruloyl-CoA. At the same time, product inhibition leads to a backward accumulation in upstream metabolites, which compensates, at least partially, for the initial decrease in p-coumaroyl-CoA. Product inhibition is easily incorporated into the GMA model (see “Methods”). Thus, in a new round of simulations, a new set of 100,000 randomly sampled parameter values was generated as before, this time accounting for product inhibition. Again, the configurations satisfying the experimental results were recorded.

Fig. 4

Substrate competition for a shared enzyme, combined with product inhibition. The accumulation of H lignin in the 4CL transgenic line calls for a regulatory mechanism that guides the flow towards the upper branch of the pathway. Direct activation or an inhibited inhibitor can achieve this result. Simulation results support the second option

Substrate competition for a shared enzyme, combined with product inhibition. The accumulation of n class="Chemical">H lignin inpan> the 4CL transgenpan>ic linpan>e calls for a regulatory mechanism that guides the flow towards the upper branch of the pathway. Direct activation or an inpan>hibited inpan>hibitor can achieve this result. Simulation results support the second option

Although the simulations showed an improvement regarding the n class="Chemical">H lignin accumulation inpan> the 4CL knpan>ockdownpan>, no topology reached the twofold inpan>crease that was reported inpan> the literature [3].

Substrate competition for shared enzymes

Several enzymes in the n class="Chemical">lignin pathway catalyze multiple reactions with slightly differenpan>t substrates, and it is reasonable to assume substrate competition for an enpan>zyme among the multiple substrates. This competition can play an important role inpan> alterinpan>g the flow of mass inpan> a mutant plant.

We explored the consequences of substrate competition with respect to the pertinent enzyme CCR. The analysis yielded the following result. If CCR favors n class="Chemical">p-coumaroyl-CoA over pan> class="Chemical">feruloyl-CoA, due to substrate competition, the flux towards H lignin is increased. In fact, simulation analysis shows that the increase in H lignin is strong enough to match the experimental data.

It could be possible that substrate competition alone would be sufficient for increased n class="Chemical">H lignin production. We tested this conjecture with a correspondinpan>g simulation, which revealed that only the combinpan>ed model with product inpan>hibition and substrate competition matches the experimenpan>tal observations. The strenpan>gth of inpan>hibition is a priori unknpan>ownpan>, but simply becomes a parameter value inpan> the GMA model (see “Methods” section). For inpan>stance, consider the pathway inpan> Fig. 4, where X2 and X6 share the same enpan>zyme for fluxes V2 and V4. Blue arrows represenpan>t the competition betweenpan> the substrates, while red arrows represenpan>t product inpan>hibition. Inpan> this case the equation for V2 becomeswhere Y2 is the enpan>zyme catalyzinpan>g the reaction (CCR).

Inhibition of 4CL in COMT knockdown transgenics

Although product inhibition and substrate competition improve the consistency between the experimental data and numerical results in CCR1 transgenic plants, the model does not match COMT knockdown data sufficiently well. Specifically, the model does not capture the observed 30 % increase in n class="Chemical">ferulic acid inpan> COMT knpan>ockdownpan>s [4]. This observationpan> becomes evenpan> more difficult to explainpan> if onpan>e conpan>siders the simultanpan>eous 20 % decrease inpan> pan> class="Chemical">coniferyl aldehyde. One could speculate that the high accumulation in 5-OH-ferulic acid might trigger a cascade of product inhibition that leads to the accumulation of ferulic acid, but computational results did not support the idea.

Further analysis with the model revealed that the reaction from n class="Chemical">ferulic acid to pan> class="Chemical">feruloyl-CoA, which is catalyzed by 4CL, is the bottleneck. Indeed, the computational results show that this reaction has a flux that is 10 times as large as the efflux from ferulic acid towards 5-OH-ferulic acid. Thus, if the flux towards ferulic acid decreases, any substantial accumulation is impossible unless the 4CL reaction is inhibited. This model-based deduction is indirectly supported by experimental data from one of our labs that exhibit a slight accumulation in the distant p-coumaric acid and caffeic acid, which is explained by 4CL inhibition as well (data not shown).

Accounting for the deduced 4CL inhibition in the model leads to simulations that faithfully capture all experimental data associated with the COMT knockdown; in particular, the 4CL substrates accumulate and the concentration of n class="Chemical">coniferaldehyde decreases, as observed. From a biochemical poinpan>t of view, onpan>e might be inpan>terested inpan> idenpan>tifyinpan>g the inpan>hibitinpan>g agenpan>t. As it was menpan>tionpan>ed earlier, the pan> class="Chemical">5-OH-ferulic acid concentration increases by 70 % in COMT knockdown plants. While the metabolite has not been identified as a substrate for 4CL, it might be reasonable to assume that it binds to 4CL in high concentrations, due to its molecular similarity, and thereby inhibit the enzyme competitively (Fig. 5). While this hypothesis remains to be experimentally validated, the same type of substrate competition with respect to 4CL has recently been proposed by others [18]. To implement 4CL inhibition in the model in the most generic manner, we simply lowered the corresponding rate constants.

Fig. 5

Parallel reactions catalyzed by 4CL. The observed simultaneous accumulation of 4CL substrates and decrease in coniferaldehyde in COMT transgenic lines can be explained with the assumption of an inhibitory effect on the reactions catalyzed by 4CL. 5-OH-ferulic acid could be a candidate for this role. Although 5-OH-ferulic acid is not a substrate for 4CL in switchgrass, it has a similar molecular shape as ferulic acid, so that high concentrations of 5-OH-ferulic acid might exert competitive inhibition that is comparable to the inhibitory effects of ferulic acid

Parallel reactions catalyzed by 4CL. The observed simultaneous accumulation of 4CL substrates and decrease in n class="Chemical">coniferaldehyde inpan> COMT tranpan>sgenpan>ic linpan>es canpan> be explainpan>ed with the assumptionpan> of anpan> inpan>hibitory effect onpan> the reactionpan>s catalyzed by 4CL. pan> class="Chemical">5-OH-ferulic acid could be a candidate for this role. Although 5-OH-ferulic acid is not a substrate for 4CL in switchgrass, it has a similar molecular shape as ferulic acid, so that high concentrations of 5-OH-ferulic acid might exert competitive inhibition that is comparable to the inhibitory effects of ferulic acid

Compatible configurations

The mathematical model with universal product inhibition, substrate competition for CCR1, inhibition of 4CL, and the possibility of a metabolic channel was subjected to large-scale simulations aimed at inferring the most likely topology of the n class="Chemical">lignin pathway (recall Fig. 3). Similar to previous simulationpan>s, a sample of 100,000 parameter sets was genpan>erated to test model conpan>sistenpan>cy with the experimenpan>tal data anpan>d to provide likely kinpan>etic orders for the model (see “Methods”). Inpan>triguinpan>gly, the onpan>ly pathway conpan>figurationpan> that is compatible with all available data is Conpan>figurationpan> 1 of Fig. 3. Note that the speculated pan> class="Chemical">coniferaldehyde channel is indeed present. In fact, no parameter set, using Configurations 2 and 3, could reproduce the experimental data which eliminates the chance to compare the relative performance of the configurations.

Principal component analysis

To gain a better understanding of the parameter space of the system, principal component analysis (PCA) was performed on the parameter sets that had been filtered by the model criteria. Once the principal components of the parameter space were identified, a new round of simulations was executed. Specifically, a sample of 100,000 parameter sets was generated along the principal directions and within the reduced space. The set was then transformed back to the original coordinates. The successful parameter sets were recorded and are depicted in Additional file 2: Figure S8. Ultimately, principal components 1 through 4 collectively account for 88 % of the variance.

Model uniqueness

It is theoretically impossible to proof the uniqueness of a model for such complex nonlinear problem, because it is always possible to evoke additional processes in such a fashion that the original model could be subsumed as a simpler special case. In our case, one should note that our large-scale simulation approach led to a structurally and numerically compact ensemble of similar solutions within the high-dimensional parameter space of the system. Given that we determined the ensemble with Monte Carlo simulations that cast a very wide net over the parameter space, it is difficult to imagine entirely different parameterizations that would capture all data as well as our ensemble and perform well in the validation studies we performed.

Moreover, considering that the available data were obtained from several independent transgenics, and that the stoichiometric system of the system is underdetermined, the likelihood of significantly other solutions appears to be rather small. Also, our simulations show that the system converges to the same steady-state starting from a wide array of initial conditions. Some arbitrary initial conditions actually lead to steady-state values outside of the defined physiological bounds; however, among the initial conditions that lead to admissible steady-states, several rounds of screening showed identical results.

In summary, it is well understood that model design is an iterative procedure, and while our logical analysis of numerical results suggested the sten class="Chemical">p-wise addition or eliminpan>ation of new features, there is no mathematical proof that the model enpan>semble is truly unique.

Outside these purely mathematical arguments, we might also look at the biological reasonableness of the model. For instance, one could ask why only CCR was subjected to substrate competition, while there are other shared enzymes. The answer is a matter of simplicity, as suggested by Ockham’s razor. Namely, we demonstrate that the substrate competition of CCR is needed to match the available data, while additional mechanisms are not necessary to explain the experimental data. Thus, we cannot exclude that additional regulatory mechanisms might exist, but we would need additional, independent data to confirm or refute such a hypothesis.

We also note that, although the model design progressed iteratively, we carefully investigated the necessity of including each individual mechanism a posteriori. For example, upon discovering that competitive inhibition over CCR improves n class="Chemical">H-lignin accumulationpan>, we asked whether product inpan>hibitionpan> was still vital for the model to explainpan> the observationpan>s. We examinpan>ed this hypothesis anpan>d determinpan>ed that pan> class="Chemical">H-lignin accumulation could not be captured anymore. We therefore concluded that both mechanisms, product inhibition and CCR competition, are necessary. We found this conclusion reasonable, as both product inhibition and substrate competition are common in metabolic pathway systems.

Model validation

The model with parameter values described above was constructed based on experimental data from wild-type n class="Species">switchgrass anpan>d four tranpan>sgenpan>ic linpan>es (4CL, CCR1, CAD anpan>d COMT knpan>ock-downpan>s). To validate the model, experimenpan>tal data from a separate tranpan>sgenpan>ic planpan>t, which had not beenpan> used inpan> anpan>y way durinpan>g the model design, were used to inpan>vestigate how well the system performs unpan>der unpan>tested conpan>ditionpan>s. Namely, inpan> a recenpan>t study, the tranpan>scriptionpan> inpan>hibitor pan> class="Chemical">PvMYB4 was over-expressed in order to reduce enzyme expression in the lignin pathway [19]. While metabolite concentrations were not measured for any of the pathway intermediates, the published data contain H, G and S lignin levels, as well as comparisons of enzyme activities between the wild type and PvMYB4 plants. The overall result of the study is a global reduction in the expression of the enzymes of the pathway, which in turn leads to 40–70 % decreases in total lignin.

We tested our model against the profile of observed enzyme expression under overexpression of n class="Chemical">PvMYB4. We started with the already parameterized model without inpan>troducinpan>g anpan>y alterationpan>s or adjustmenpan>ts, except for resettinpan>g the appropriate enpan>zyme activities, anpan>d tested how the system responpan>ded to the inpan>hibitionpan> inpan> comparisonpan> to the inpan> vivo experimenpan>ts [19]. Enpan>couraginpan>gly, the altered G- anpan>d S-pan> class="Chemical">lignin amounts and their ratio, reported in the experimental study, are captured by the model with the compatible topological configuration quite well. The H-lignin was essentially unchanged in the experiment, while it slightly decreases in our model, in accordance with the data we used. However, H-lignin constitutes only about 3 % of the total lignin so that this difference is of no particular pertinence. Results are shown in Figs. 6 and 7. Figure 6 compares the fold change in lignin monomers between the experimental data and model results. The first row shows the fold change in G, S, the total lignin, and the S/G ratio comparing the wild type and PvMYB4 lines from the experiment; the second row corresponds to the computed configuration. As can be seen, the model results are quite consistent with experimental data.

Fig. 6

Fold changes in lignin monomer concentrations in PvMYB4 transgenic plants. The top row represents the average of PvMYB4 plants experimental data normalized with respect to the average of the control plants. The second row represent the results of the model with settings corresponding to the PvMYB4 experiment in [19], normalized with respect to wild-type model results. Wild type is set to 1, which corresponds to white in the color bar. H lignin only counts for 3 % of total lignin and is not shown in here

Fig. 7

Steady-state profiles of key pathway metabolites in PvMYB4 overexpression as predicted by the model. Concentrations are normalized and the base value is set to 100, which corresponds to white in the color bar. Any increases with respect to the wild-type steady state are reflected in the red spectrum and any decreases in the blue spectrum

Fold changes in n class="Chemical">lignin monpan>omer conpan>cenpan>trationpan>s inpan> pan> class="Chemical">PvMYB4 transgenic plants. The top row represents the average of PvMYB4 plants experimental data normalized with respect to the average of the control plants. The second row represent the results of the model with settings corresponding to the PvMYB4 experiment in [19], normalized with respect to wild-type model results. Wild type is set to 1, which corresponds to white in the color bar. H lignin only counts for 3 % of total lignin and is not shown in here

Steady-state profiles of key pathway metabolites in n class="Chemical">PvMYB4 overexpression as predicted by the model. Concenpan>trations are normalized and the base value is set to 100, which corresponds to white inpan> the color bar. Any inpan>creases with respect to the wild-type steady state are reflected inpan> the red spectrum and any decreases inpan> the blue spectrum

This independent validation is very reassuring, especially with respect to future attempts to use metabolic engineering techniques to alter the S/G ratio in n class="Species">switchgrass. For inpan>stance, if further model predictions prove similarly reliable, the model could be used to simulate and optimize the outcome of combinpan>atorial knpan>ockdownpan>s, whose outcomes are not necessarily predictable with inpan>tuition alone. Such predictions would be very valuable, as a comprehenpan>sive combinpan>atorial screenpan>inpan>g of double and triple knpan>ock-downpan>s would neither be economical nor experimenpan>tally feasible.

While the published PvMYB4 data used for the first validation do not containpan> inpan>termediate metabolite concenpan>trations, a more recenpan>t study provides steady-state data for several of the pathway metabolites [20]. Comparinpan>g the published data inpan> [20] with those inpan> our model, we finpan>d that sevenpan> metabolites are represenpan>ted inpan> both, namely, caffeic acid, 5-OH-coniferyl alcohol, ferulic acid, sinapyl alcohol, coniferaldehyde, p-coumaric acid and coniferyl alcohol.

Figure 7 exhibits a comparison of the steady-state profiles. The top row shows the simulation results, while the bottom row represents experimentally measured steady-state concentrations in n class="Chemical">PvMYB4 normalized to wild type from [20]. The wild-type value for each conpan>cenpan>trationpan> is set to 100 (white), anpan>d the red-blue spectrum represenpan>ts inpan>creases or decreases inpan> steady-state values of knpan>ockdownpan>s. For five of these sevenpan> metabolites, our computationpan>al results of pan> class="Chemical">PvMYB4 conditions show the same semi-quantitative behavior in steady-state concentrations compared to the wild type; these are caffeic acid, 5-OH-coniferyl alcohol, sinapyl alcohol, coniferaldehyde and coniferyl alcohol. Discrepancies are seen in ferulic acid and p-coumaric acid. Here, the experimental data show a decrease in the steady-state concentrations, while our computational results predict an accumulation. Interestingly, these differences occur for metabolites whose effluxes out of the lignin pathway are ill defined, because their characteristics were not documented in the literature. It is therefore likely that they are not optimally parameterized in the model.

Discussion

In this work, we developed an ensemble of models of n class="Chemical">lignin biosynpan>thesis inpan> stem anpan>d tiller tissue inpan> pan> class="Species">switchgrass, P. virgatum. The model reflects the consequences of various enzyme knock-downs quite well and performed satisfactorily in two validation studies with experimental data that had not been used in the model design or implementation. We used as the modeling framework the generalized mass action (GMA) format within biochemical systems theory (BST) [21-25]. The power-Law representation, which is the hallmark of this type of model, is arguably the least biased default formulation and by its mathematical nature avoids problems due to possibly invalid assumptions that may cast doubt on traditional Michaelis–Menten models in vivo [26]. Parameter values were, as always, difficult to obtain in a direct manner. We used for this purpose experimental knock-down data and a sophisticated Monte Carlo sampling strategy that has been used very successfully for similar systems before [14]. As a particular sub-goal, we investigated the regulatory mechanism of the pathway and the possible co-localization or coupling of the pair of enzymes, CCR1/CAD that was previously suggested for Medicago [5].

To elucidate the co-localization or coupling of these enzymes in n class="Species">switchgrass, we studied multiple conpan>figurationpan>s that seemed a priori plausible anpan>d idenpan>tified those natural designs that were conpan>sistenpan>t with the experimenpan>tal data. The conpan>sistenpan>t designs were further examinpan>ed unpan>der differenpan>t regulationpan> scenpan>arios. The mainpan> result from this study is a very robust model of pan> class="Chemical">lignin biosynthesis in switchgrass that is consistent with all available data. The model was, at least to some degree, validated with a formerly unused dataset. If this validation can be confirmed and expanded experimentally, the model proposed here may be used to predict responses of the natural pathway system to alterations that are difficult to assess with experimental means. For instance, a further validated model will allow the prediction of responses to combinatorial knockdowns that could be the basis for future designs of more sophisticated transgenic lines than are currently available.

The computational analysis suggests the co-localization or functional coupling of the two enzymes CCR1 and CAD. Metabolic channeling and compartmentalization in plants have been identified in many biochemical pathways [27]. Of importance here, it has been suggested that enzymes catalyzing early reactions in the n class="Chemical">monolignol pathway may be co-localized inpan> their binpan>dinpan>g to the ER. For inpan>stanpan>ce, a multi-proteinpan> complex has beenpan> idenpan>tified betweenpan> PAL anpan>d pan> class="Chemical">C4H, and it seems that most of the substrates use these channels, but that some substrate undergoes the metabolic conversion in two steps [28-30]. C4H can also form a complex with C3′H [31], and it has been suggested that different forms of 4CL form a complex in poplar [32]. Independent computational work on alfalfa came to a similar conclusion for channeling of enzymes associated with coniferaldehyde, which were proposed to form a metabolic channel [14]. Our results on switchgrass, presented in this article, are in line with the latter result and suggest moreover that channeling around coniferaldehyde is necessary to capture the available data.

The comparative study of different configurations revealed that consistency with the available experimental data was most difficult to achieve for transgenic 4CL down-regulated lines, in which, surprisingly, the n class="Chemical">H lignin conpan>cenpan>trationpan> is inpan>creased. This observationpan> is at first counpan>terinpan>tuitive because 4CL is located directly upstream of the pan> class="Chemical">H lignin precursors, which would lead to the a priori expectation of a decrease in H lignin. The combination of two postulated types of regulatory mechanisms was able to explain this observation. The first is product inhibition, which is observed quite frequently in biochemical systems. While improving the data compatibility, this mechanism turned out to be insufficient, thus requiring additional signaling. Arguably the simplest explanation is a regulatory structure that works in either of the mechanisms below:

An intermediate in the pathway is increased in response to the 4CL knockdown and activates the precursors of H lignin synthesis. The most likely candidates for this scenpan>ario appear to be p-coumaric acid, caffeic acid, and ferulic acid (Fig. 8a).

Fig. 8

Two plausible explanations for an increase in the H lignin concentration in 4CL transgenic lines. a represents a putative increase in an activator located upstream of the enzyme 4CL, whereas b shows a putative decrease in an inhibitor located downstream of 4CL

Two plausible explanations for an increase in the n class="Chemical">H lignin concenpan>tration inpan> 4CL transgenpan>ic linpan>es. a represenpan>ts a putative inpan>crease inpan> an activator located upstream of the enpan>zyme 4CL, whereas b shows a putative decrease inpan> an inpan>hibitor located downpan>stream of 4CL

There exists an inhibitor for the H lignin branch. This metabolite would have to be located such that its concenpan>tration is decreased due to the 4CL knpan>ockdownpan>, which means that the inpan>hibitor activity is inpan>hibited and therefore exerts a net positive effect on the system (Fig. 8b). Feruloyl-CoA could be a good candidate for this scenario.

The current literature does not support the first hypothesis. By contrast, multiple candidates are available for the second scenario. A reasonable scenario arises from the fact that the n class="Chemical">lignin pathway inpan> pan> class="Species">switchgrass includes parallel fluxes that share the same enzymes. Indeed, 4CL, CAD, COMT, F5H and CCR1 all catalyze multiple reactions, and it is likely that the substrates exert competitive inhibition for the shared enzyme, as it was also suggested in [33]. Supporting this scenario, a targeted numerical analysis demonstrated that competition over CCR1 perfectly matches the results of the 4CL knockdown line in the model with product inhibition. One could surmise that the latter mechanism would suffice to represent the increase in H lignin concentration. To test this hypothesis, we simulated the model with enzyme competition but without product inhibition. The results showed that competitive inhibition by itself could not satisfactorily resolve the issue. By contrast, the combined model containing product inhibition and competitive inhibition matches the experimental results very well. One should also recall that the product inhibition and substrate competition mechanisms only work properly if the proposed metabolic channel is present (Fig. 3, Configuration 1).

Another aspect of the experimental data that was not captured well by the original model, even when product inhibition and substrate competition over CCR1 were taken into account, is the accumulation of 4CL substrates in COMT transgenic plants. Particularly counterintuitive appears to be the accumulation of n class="Chemical">ferulic acid as a product of a reactionpan> catalyzed by COMT. The observed conpan>comitanpan>t decrease inpan> the steady-state conpan>cenpan>trationpan> of pan> class="Chemical">coniferaldehyde supports the possible explanation that the observation is due to regulation that begins to inhibit the conversion of ferulic acid into coniferaldehyde, when 4CL substrates are in excess. The simultaneous accumulation of p-coumaric acid and caffeic acid provides additional evidence that reactions catalyzed by 4CL are inhibited in COMT knockdown plants. Accounting for this feature to our model, all experimental data are represented well. The mechanism of the regulation remains a subject of further experimental investigations. Figure 9 shows the pathway including all inferred regulatory signals.

Fig. 9

Full scheme of the lignin biosynthetic pathway in switchgrass suggested by the computational results of this study. All regulatory signals, i.e., universal product inhibition, substrate competition over CCR1, and 4CL inhibition are shown. The 4CL inhibiting agent is unknown and therefore denoted with X. 5-OH-ferulic acid might be a candidate for this role

Full scheme of the n class="Chemical">lignin biosynpan>thetic pathway inpan> pan> class="Species">switchgrass suggested by the computational results of this study. All regulatory signals, i.e., universal product inhibition, substrate competition over CCR1, and 4CL inhibition are shown. The 4CL inhibiting agent is unknown and therefore denoted with X. 5-OH-ferulic acid might be a candidate for this role

Conclusions

The model proposed in this article captures all available data and performed well in independent n class="Chemical">PvMYB4 validationpan> experimenpan>ts. This good match with data is reasonpan> for cautious optimism, which however is to be supported with further experimenpan>tal conpan>firmationpan>. Inpan>deed, work is inpan> progress to genpan>erate anpan>d anpan>alyze additionpan>al tranpan>sgenpan>ic pan> class="Species">switchgrass lines and to incorporate further lignin compositional and enzyme activity and kinetic data into the model. If the model fares well in these additional validation studies, the results from the present study suggest that one might use the model for predictions, for instance, with respect to double knock-downs, and for optimization studies that could potentially affect the lignin-based recalcitrance in switchgrass in a favorable manner.

Methods

Model construction

Much of the analysis in this article consists of comparisons and simulations with different models. Each of these models consists of a system of differential equations that represent the rate of change in metabolite concentrations, which are represented as dependent variables. The right-hand side of each equation contains a set of fluxes which enter (influxes) or leave (effluxes) the metabolite pool. Enzymes are included in the model as independent variables; that is, they do not change in activity during any given computational experiment. The generic formulation of each equation iswhere each X is a metabolite, V are fluxes associated with X, and the quantities s are stoichiometric coefficients, which here are simply 0, 1 or −1 and determine whether flux V affects X as influx or efflux or not at all. Each V is a function of some or potentially all of the X. At the steady state, the left-hand side is equal to zero, and fluxes can be assessed with methods of linear algebra [34]. Because the system in our case is underdetermined, infinitely many solutions satisfy the steady-state condition. Following the tenets of Flux Balance Analysis (FBA), an objective function is chosen and the problem is solved as a linear programming problem [34]. In the present study, maximizing the total amount of n class="Chemical">lignin is set as the objective of the system. The optimizationpan> problem is solved usinpan>g MATLAB (versionpan> R2014a, The MathWorks, Natick, MA, USA) funpan>ctionpan> linpan>pro. The output is the set of fluxes at the steady state that maximizes the definpan>ed objective.

The fluxes themselves are formulated as general mass action (GMA) models of the typewithin the modeling framework of BST [21, 22, 24, 35, 36]. Here, α is the rate constant, each X, for 1 < r < n, is a metabolite or, for n + 1 < r < n + m, an enzyme involved in the reaction. Thus, n is the number of metabolites and m is the number of enzymes in the pathway. The exponents g are kinetic orders that quantify the effect of X on V. Similarly, h describes the effect of the enzyme on the reaction. It is customary to set each h to 0 or 1, thus merely reflecting absence or presence of an enzyme in a specific flux. This setting of h = 1 is consistent with the underpinnings of Michaelis–Menten, mass-action, and other traditional models, where a reaction is assumed to be a linear function of enzyme activity. All other kinetic orders g are sampled from the range between 0 and 1 if X is a substrate or activator of the flux, or from the range between −1 and 0 if X is an inhibitor.

Due to the nature of the present experimental data for n class="Species">switchgrass, the real concenpan>trations of metabolites and enpan>zyme activities inpan> vivo are unknpan>ownpan>. As a remedy, we normalize these quantities with respect to the steady state and set all base values to 100. Thus, we setand express Eq. (4) as

Since the constant X refers to the steady state, simple algebra adjusts the rate constants to this steady state. Thus, we obtain

The enzymes are independent variables and therefore constant for each experiment. Therefore, X = X for n + 1 < r < n + m for wild type, whereas for a transgenic line it takes a value between 0 and 1, according to the level of knockdown. At the steady state we have:with this setting, each steady-state flux is given as

If the flux is known, the rate constant can be computed as

With these settings, the set of the differential equations for the model takes the form below.

where the quantities include the influxes into the pathway and the fluxes, , are defined as follows:

The metabolites of the pathway arewhile the enzymes of the pathway are

Note that the model does not account for the dynamics of n class="Chemical">tyrosine, which we consider constant here. The model scheme is shownpan> inpan> Fig. 10.

Fig. 10

Lignin pathway in the notation of the model. Redundancy of enzymes, i.e., 4CL, CCR1, CAD, COMT and F5H in parallel fluxes reduces the dimension of state space. The enzymes HCT, C3′H and CSE in flux V9 are merged into one independent variable, Z25. Note that the presence of the G-channel, V26, is an inference from the computational simulations results

n class="Chemical">Lignin pathway inpan> the notationpan> of the model. Redunpan>danpan>cy of enpan>zymes, i.e., 4CL, CCR1, CAD, COMT anpan>d pan> class="Chemical">F5H in parallel fluxes reduces the dimension of state space. The enzymes HCT, C3′H and CSE in flux V9 are merged into one independent variable, Z25. Note that the presence of the G-channel, V26, is an inference from the computational simulations results

Parameter space and sampling

Similar to earlier work [14, 37, 38], flux rates are computed with FBA. Next, the parameters to be estimated are the kinetic orders and rate constants are in turn estimated from the FBA results and randomly sampled kinetic orders through the steps mentioned above. The kinetic order of a metabolite is positive if the metabolite is a substrate or activator of the flux and negative if it acts as an inhibitor. The kinetic order of each enzyme has a default value of 1, which is in line with traditional enzyme kinetics, because it is customary to assume that a flux has a linear relationship with the enzyme. This assumption is explicitly or implicitly made in essentially all traditional models of enzyme kinetics as, for instance, in the Michaelis–Menten formalism, where Vmax equals kcat times the enzyme concentration.

The down-regulation of an enzyme is modeled through the enzyme concentration, not the kinetic order. Since the concentrations of metabolites and enzymes are normalized, the concentration of an enzyme in the wild type has the default value of 1. In transgenics, the concentration of the corresponding enzyme is set to a value less than one if it is down-regulated. For example, to represent the 4CL knockdown, the concentration of the enzyme is set to 0.6 as the enzyme is down-regulated by 40 %.

To account for product inhibition, the inhibiting product is represented in each reaction by a factor consisting of its concentration, raised to a negative power. The result is as follows:

Here, S is the substrate, P is the product, g is the kinetic order of the inhibiting product and g is the kinetic order of the substrate. The ratio of kinetic orders could be derived directly [22] from the corresponding expression for a product-inhibited Michaelis–Menten reaction, which takes the form

The power-law form of Eq. 15 can directly be computed from the tenets of Biological Systems Theory (BST), which defines the kinetic orders as

Rearrangement of these equations gives the ratio of kinetic orders as follows:

The bounded ratio of kinetic orders provides a valuable constraint for the Monte Carlo simulations, because a fixed ratio does not affect the dimension of the parameter space.

For the initial set of simulations, the sampling space is chosen as a unit hypercube in ℝ where n is number of kinetic orders to be estimated. A set of 100,000 parameter sets is generated for each scenario simulation. 10,000 sets are randomly generated from the sampling space using Latin Hypercube Sampling to assure a homogeneous coverage of the space, while 90,000 sets are generated by the MATLAB (version R2014a, The MathWorks, Natick, MA, USA) function rand. Each parameter set is simulated to examine whether the model with this set can match the experimental results for the wild type and transgenics. The model is deemed a match for the experimental results if:

The model returns proper n class="Chemical">lignin contenpan>ts and S/G ratios for the wild type and differenpan>t transgenpan>ics, with downpan>-regulation of 4CL (40 %), CCR1 (50 %), COMT (30 %), and CAD (30 %).

The model returns the proper decrease in n class="Chemical">lignin contenpan>t inpan> the case of knpan>ockdownpan>s inpan> 4CL, CCR1, COMT, and CAD.

The model demonstrates an increase in n class="Chemical">H lignin inpan> 4CL transgenpan>ics.

The model matches the altered metabolite concentrations in the COMT transgenic.

If a parameter satisfies the above conditions, it is recorded along with the corresponding topological configuration.

While our model approach emphasizes ensembles of feasible models, the parameter values in Tables 2, 3, and 4 represent one implementation, which we used for further numerical exploration. This specific parameter set corresponds to the minimum error in the comparison of the model results in n class="Chemical">PvMYB4 anpan>d the experimenpan>tal data.

Table 2

A sample of rate constants from the ensemble of rate constants

α_1	0.5233	α_8	0.0058	α_15	0.0771	α_22	0.0392
α_2	0.1053	α_9	0.2265	α_16	0.0881	α_23	0.1573
α_3	0.15	α_10	0.0024	α_17	0.0168	α_24	0.0712
α_4	0.2711	α_11	0.1054	α_18	0.002	α_25	0.1154
α_5	0.1832	α_12	0.1095	α_19	0.2212	α_26	0.0814
α_6	0.003	α_13	0.1452	α_20	0.0402
α_7	0.0042	α_14	0.2681	α_21	0.0002

Table 3

A sample of kinetic orders from the ensemble of kinetic orders

g_1,1	0.2813	g_6,7	0.4040	g_10,14	−0.1023	g_14,21	0.8535
g_2,1	−0.1406	g_3,8	0.5759	g_10,15	0.9009	g_13,22	0.7673
g_2,2	0.0846	g_4,9	0.6118	g_11,15	−0.4505	g_15,22	−0.3836
g_2,3	0.8240	g_7,9	−0.3059	g_4,15	−0.0355	g_14,23	0.1043
g_3,3	−0.4120	g_4,10	0.6398	g_11,16	0.5198	g_16,23	−0.0521
g_3,4	0.5669	g_7,11	0.6277	g_12,16	−0.2599	g_15,24	0.5080
g_4,4	−0.2835	g_8,11	−0.3138	g_12,17	0.7121	g_16,24	−0.2540
g_4,5	0.0710	g_7,12	0.6414	g_9,18	0.6982	g_16,25	0.2855
g_5,5	−0.0355	g_9,12	−0.3207	g_11,19	0.0160	g_10,26	0.7116
g_10,5	−0.4505	g_8,13	0.4885	g_13,19	−0.0080	g_12,26	−0.3558
g_5,6	0.9630	g_10,13	−0.2442	g_12,20	0.6973	g_4,26	−0.0355
g_6,6	−0.4815	g_9,14	0.2046	g_14,20	−0.3487

Table 4

Initial values

Z_0,1	100	Z_0,5	100	Z_0,9	100	Z_0,13	100
Z_0,2	100	Z_0,6	100	Z_0,10	100	Z_0,14	100
Z_0,3	100	Z_0,7	100	Z_0,11	100	Z_0,15	100
Z_0,4	100	Z_0,8	100	Z_0,12	100	Z_0,16	100