Glycan processing in the Golgi -- optimal information coding and constraints on cisternal number and enzyme specificity

Many proteins that undergo sequential enzymatic modification in the Golgi cisternae are displayed at the plasma membrane as cell identity markers. The modified proteins, called glycans, represent a molecular code. The fidelity of this glycan code is measured by how accurately the glycan synthesis machinery realises the desired target glycan distribution for a particular cell type and niche. In this paper, we quantitatively analyse the tradeoffs between the number of cisternae and the number and specificity of enzymes, in order to synthesize a prescribed target glycan distribution of a certain complexity. We find that to synthesize complex distributions, such as those observed in real cells, one needs to have multiple cisternae and precise enzyme partitioning in the Golgi. Additionally, for fixed number of enzymes and cisternae, there is an optimal level of specificity of enzymes that achieves the target distribution with high fidelity. Our results show how the complexity of the target glycan distribution places functional constraints on the Golgi cisternal number and enzyme specificity.


I. INTRODUCTION
A majority of the proteins synthesized in the endoplasmic reticulum (ER) are transferred to the Golgi cisternae for further chemical modification by glycosylation [1], a process that sequentially and covalently attaches sugar moieties to proteins, catalyzed by a set of enzymatic reactions within the ER and the Golgi cisternae.These enzymes, called glycosyltransferases, are localized in the ER and cis-medial and trans Golgi cisternae in a specific manner [2,3].Glycans, the final products of this glycosylation assembly line are delivered to the plasma membrane (PM) conjugated with proteins, whereupon they engage in multiple cellular functions, including immune recognition, cell identity markers, cell-cell adhesion and cell signalling [2][3][4][5][6].This glycan code [7,8], representing information [9] about the cell, is generated dynamically, following the biochemistry of sequential enzymatic reactions and the biophysics of secretory transport [4,10,11].
In this paper, we will focus on the role of glycans as markers of cell identity.For the glycans to play this role, they must inevitably represent a molecular code [4,7,11].While the functional consequences of glycan alterations have been well studied, the glycan code has remained an enigma [7,[11][12][13].In this paper, we study one aspect of molecular coding, namely the fidelity of this molecular code generation.While it has been recognised that fidelity of the glycan code is necessary for reliable cellular recognition [14], a quantitative measure of fidelity of the code and its implications for cellular structure and organization are lacking.
There are two aspects of the cell-type specific glycan code that have an important bearing on quantifying fidelity.The first is that extant glycan distributions have high complexity, owing to evolutionary pressures arising from (a) reliable cell type identification amongst a large set of different cell types in a complex organism, the preservation and diversification of "self-recognition" [5], (b) pathogen-mediated selection pressures [2,4,6], and (c) herd immunity within a heterogenous population of cells of a community [15] or within a single organism [5].Here, we will interpret this to mean that the target distribution of glycans of a given cell type is complex; in Sect.II we define a quantitative measure for complexity and demonstrate its implications in the context of human T-cells.The second is that the cellular machinery for the synthesis of glycans, which involves sequential chemical processing via cisternal resident enzymes and cisternal transport, is subject to variation and noise [4,10,11]; the synthesized glycan distribution is, therefore, a function of cellular parameters such as the number and specificity of enzymes, inter-cisternal transfer rates, and number of cisternae.We will discuss an explicit model of the cellular synthesis machinery in Sect.III.
In this paper, we define fidelity as the minimum achievable Kullback-Leibler (KL) divergence [16,17] between the synthesized distribution of glycans and the target glycan distribution.This KL divergence is a function of the cellular parameters governing glycan synthesis, such as the number and specificity of enzymes, inter-cisternal transfer rates, and number of cisternae (Sect.V).We analyze the tradeoffs between the number of cisternae and the number and specificity of enzymes, in order to achieve a prescribed target glycan distribution with high fidelity (Sect.VI).Our analysis leads to a number of interesting results, of which we list a few here: (i) In order to construct an accurate representation of a complex target distribution, such as those observed in real cells, one needs to have multiple cisternae and precise enzyme partitioning.Low complexity target distributions can be achieved with fewer cisternae.(ii) This definition of fidelity of the glycan code, allows us to provide a quantitative argument for the evolutionary requirement of multiple-compartments.(iii) For fixed number of enzymes and cisternae, there is an optimal level of specificity of enzymes that achieves the complex target distribution with high fidelity.(iv) Keeping the number of enzymes fixed, having low specificity or sloppy enzymes and larger cisternal number could give rise to a diverse repertoire of functional glycans, a strategy used in organisms such as plants and algae.
Stated another way, our results imply that the pressure to achieve the target glycan code for a given cell type, places strong constraints on the cisternal number and enzyme specificity [18].This would suggest that a description of the nonequilibrium assembly of a fixed number of Golgi cisternae must combine the dynamics of chemical processing with membrane dynamics involving fission, fusion and transport [19,20], opening up a new direction for future research.

II. COMPLEXITY OF GLYCAN CODE IN REAL CELLS
Since each cell type (in a niche) is identified with a distinct glycan profile [4,7,11], and this glycan profile is noisy because of the stochastic noise associated with the synthesis and transport [11][12][13], a large number of different cell types can be differentiated only if the cells are able to produce a large set of glycan profiles that are distinguishable in the presence of this noise.A more complex or richer class of glycan profiles is able to support a larger number of well separated profiles, and therefore, a larger number cell types, or equivalently, a more complex organism 1In order to implement a quantitative measure of complexity, we first need a consistent way of smoothening or coarsegraining the discrete glycan distribution to remove measurement and synthesis noise.In this paper, we approximate the glycan profile as mixture of Gaussian densities with specified number of components that are supported on a finite set of indices [21].Since the complexity of k-component Gaussian is an increasing function of k, we use the number of component k and complexity interchangeably.
Using this definition we demonstrate that the glycan profiles of typical mammalian cells are very complex.We obtain target profiles for a given cell type from Mass Spectrometry coupled with determination of molecular structure (MSMS) measurements [22].Fig. 1 shows the the MSMS data from human T-cells and human and mouse neutrophils [22], and their coarse-grained representations using Gaussian mixture models (GMM) of differing complexity -a low complexity k = 3 GMM and high complexity k = 20 GMM.It is clear from Fig. 1, that the more complex k = 20 GMM is a better representation of the MSMS data as compared to the less complex k = 3 GMM.Indeed the k = 20 Gaussian mixture model is the best compromise between faithfulness of the representation and cost of an additional component, as seen from the saturation of the likelihood function [17].Details of this systematic coarse-graining procedure appear in Sect.VI B and Appendix G.
Having demonstrated the complexity of the typical glycan distributions associated with a given cell type, we will now describe a general model of the cellular machinery that is capable of synthesizing glycans of the desired complexity.We expect that cells need a more elaborate mechanism to produce profiles from a more complex set.

III. SYNTHESIS OF GLYCANS IN THE GOLGI CISTERNAE
The glycan display at the cell surface is a result of proteins that flux through and undergo sequential chemical modification in the secretory pathway, comprising an array of Golgi cisternae situated between the ER and the PM, are injected from the ER to cisterna-1 at rate q.Superimposed is transition network of chemical reactions (column) -intercisternal transfer (rows), the latter with rates µ (j) .Pc (j) k denotes the acceptor substrate in compartment j and the glycosyl donor c0 is chemostated in each cisterna.This results in a frequency distribution of glycans displayed at the PM (red curve), that is representative of the cell type.
as depicted in Fig. 2. Glycan-binding proteins (GBPs) are delivered from the ER to the first cisterna, whereupon they are processed by the resident enzymes in a sequence of steps that constitute the N-glycosylation process [2].A generic enzymatic reaction in the cisterna involves the catalysis of a group transfer reaction in which the monosaccharide moiety of a simple sugar donor substrate, e.g.UDP-Gal, is transferred to the acceptor substrate, by a Michaelis-Menten (MM) type reaction [2] Acceptor + glycosyl donor + Enzyme From the first cisterna, the proteins with attached sugars are delivered to the second cisterna at a given inter-cisternal transfer rate, where further chemical processing catalysed by the enzymes resident in the second cisterna occurs.This chemical processing and inter-cisternal transfer continues until the last cisterna, thereupon the fully processed glycans are displayed at the PM [2].The network of chemical processing and inter-cisternal transfer forms the basis the physical model that we will describe next.
Any physical model of such a network of enzymatic reactions and cisternal transfer needs to be augmented by reaction and transfer rates and chemical abundances.To obtain the range of allowable values for the reaction rates and chemical abundances, we use the elaborate enzymatic reaction models, such as the KB2005 model [23][24][25] (with a network of 22, 871 chemical reactions and 7565 oligosaccharide structures) that predict the N-glycan distribution based on the activities and levels of processing enzymes distributed in the Golgi-cisternae of mammalian cells, and compare these predictions with N-glycan mass spectrum data.For the allowable rates of cisternal transfer, we rely on the recent study by Ungar and coworkers [26,27], whose study shows how the overall Golgi transit time and cisternal number, can be tuned to engineer a homogeneous glycan distribution.

A. Chemical reaction and transport network in cisternae
With this background, we now define our quantitative model for chemical processing and transport in the secretory pathway.We consider an array of N C Golgi cisternae, labelled by j = 1, . . ., N C , between the ER and the PM (Fig. 2).
Glycan-binding proteins (GBPs), denoted as Pc 1 , are delivered from the ER to cisterna-1 at an injection rate q.It is well established that the concentration of the glycosyl donor in the j-th cisterna is chemostated [2,[28][29][30], thus in our model we hold its concentration c (j) 0 constant in time for each j.The acceptor Pc    3 , and so on, provided the requisite enzymes are present in that cisterna.This leads to the sequence of enzymatic reactions Pc k → . .., where k enumerates the sequence of glycosylated acceptors, using a consistent scheme (such as in [23]).The glycosylated GBPs are transported from cisterna-1 to cisterna-2 at an inter-cisternal transfer rate µ (1) , whereupon similar enzymatic reactions proceed.The processes of intra-cisternal chemical reactions and inter-cisternal transfer continue to the other cisternae and form a network as depicted in Fig. 2. Although, in this paper, we focus on a sequence of reactions that form a line-graph, the methodology we propose extends to tree-like reaction sequences, and more generally to reaction sequences that form a directed acyclic graphs [31].
Let N s denote the maximum number of possible glycosylation reactions in each cisterna j, catalysed by enzymes labelled as E (j) α , with α = 1, . . ., N E , where N E is the total number of enzyme species in each cisterna.Since many substrates can compete for the substrate binding site on each enzyme, one expects in general that N s N E .The configuration space of the network Fig. 2 is N s × N C .For the N-glycosylation pathway in a typical mammalian cell, N s = 2 × 10 4 , N E = 10 − 20 and N C = 4 − 8 [23][24][25]27].We account for the fact that the enzymes have specific cisternal localisation, by setting their concentrations to zero in those cisternae where they are not present.

Now the action of enzyme E (j)
α on the substrate Pc (j) k in cisterna j is given by In general, the forward, backward and catalytic rates ω f , ω b and ω c , respectively, depend on the cisternal label j, the reaction label k, and the enzyme label α, that parametrise the MM-reactions [32].For instance, structural studies on glycosyltransferase-mediated synthesis of glycans [33], would suggest that the forward rate ω f to depend on the binding energy of the enzyme E α to acceptor substrate Pc k and a physical variable that characterises cisterna-j.
A potential candidate for such a cisternal variable is pH [34], whose value is maintained homeostatically in each cisterna [35]; changes in pH can affect the shape of an enzyme (substrate) or their charge properties, and in general the reaction efficiency of an enzyme has a pH optimum [32].Another possible candidate for a cisternal variable is membrane bilayer thickness [36] -indeed both pH [37] and membrane thickness are known to have a gradient across the Golgi cisternae.We take ω f (j, k, α) ∝ P (j) (k, α), where P (j) (k, α) ∈ (0, 1), is the binding probability of enzyme E α with substrate Pc (j) k , and define the binding probability P (j) (k, α) using a biophysical model, similar in spirit to the Monod-Wyman-Changeux model of enzyme kinetics [38,39], of enzyme-substrate induced fit.
Let l (j) α and l k denote, respectively, the optimal "shape" for enzyme E (j) α and the substrate Pc (j) k .We assume that the mismatch (or distortion) energy between the substrate k and enzyme E (j) α , with a binding probability given by, where . is a distance metric defined on the space of l (j) α (e.g., the square of the 2 -norm would be related to an elastic distortion model [40]) and the vector σ ≡ [σ (j) α ] parametrises enzyme specificity.A large value of σ (j) α indicates a highly specific enzyme, a small value of σ (j) α indicates a promiscuous or sloppy enzyme.It is recognised that the degree of enzyme specificity or sloppiness is an important determinant of glycan distribution [2,[41][42][43].
As in [23][24][25], our synthesis model is mean-field, in that we ignore stochasticity in glycan synthesis that may arise from low copy numbers of substrates and enzymes, multiple substrates competing for the same enzymes, and kinetics of inter-cisternal transfer.Then the usual MM-steady state condition on (2), which assumes that the concentration of the intermediate enzyme-substrate complex does not change with time, implies where c k is the concentration of the acceptor substrate Pc Together with the constancy of the total enzyme concentration, E , this immediately fixes the kinetics of product formation (not including inter-cisternal transport), where From the above, the experimentally measurable parameters V max and MM-constant K M , for each (j, k, α) can be easily read out.As is the usual case, the maximum velocity V max is not an intrinsic property of the enzyme, because it is dependent on the enzyme concentration E (j) α tot ; while K M (j, k, α) = M (j, k, α)c (j) 0 /P (j) (k, α) is an intrinsic parameter of the enzyme and the enzyme-substrate interaction.The enzyme catalytic efficiency, the so-called "k cat /K M " ∝ P (j) (k, α) and is high for perfect enzymes [44] with minimum mismatch.
We now add to this chemical reaction kinetics, the rates of injection (q) and inter-cisternal transport µ (j) from the cisterna j to j + 1; in Appendix A we display the complete set of equations that describe the changes in the substrate concentrations c (j) k with time.These kinetic equations automatically obey the conservation law for the protein concentration (p).Rescaling the kinetic parameters in terms of the injection rate q, i.e.V (j, k, α) = V (j, k, α)/q and µ (j) = µ (j) /q, we see that the steady state concentrations of the glycans in each cisterna satisfy the following recursion relations (see, Appendix A).In the first cisterna, µ (1)  and in cisternae j ≥ 2, Equations ( 5)-( 6) automatically imply that the total steady state glycan concentration in each cisterna j = 1, . . ., N c is given by These nonlinear recursion equations ( 5)-( 6) have to be solved numerically to obtain the steady state glycan concentrations, c ≡ c the transport rates µ ≡ [µ (j) ] and specificity, σ ≡ [σ Now, with both the protocol for determining the target glycan distribution and the sequential chemical processing model in hand, we can precisely define the optimization problem referred to in the Introduction.Let c * denote the "target" concentration distribution2 for a particular cell type, i.e. the goal of the sequential synthesis mechanism described in Sect.IV A is to approximate c * .Let c denote the steady state glycan concentration distribution displayed on the PM - (6) We measure the fidelity between the c * and c by the Kullback-Leibler metric [16,17], Thus, the problem of designing a sequential synthesis mechanism that approximates c * for a given enzyme specificity σ, transport rate µ, number of enzymes N E , and number of cisternae N C is given by There is separation of time scales implicit in Optimization A -the chemical kinetics of the production of glycans and their display on the PM happens over cellular time scales, while the issues of tradeoffs and changes of parameters are driven over evolutionary timescales.
Optimization A, though well-defined, is a hard problem, since the steady state concentrations (6) are not explicitly known in terms of the parameters (M, V, L).In Appendix B, we formulate an alternative problem Optimization B in which the steady state concentrations are defined explicitly in terms of a new parameters R and L, and in Appendix C we prove that Optimization A and Optimization B are exactly equivalent.This is a crucial insight that allows us to obtain all the results that follow.
In Appendix H, we describe the variant of the Sequential Quadratic Programming (SQP) [45], that we use to numerically solve the optimization problem.

VI. RESULTS OF OPTIMIZATION
To start with, the dimension of the optimization search space is extremely large To make the optimization search more manageable, we ignore the k-dependence of the vectors (M, V), (or, alternatively of R, see Appendix B for details).The dependence on the reaction rates on the glycosyl substrate is still present in the forward reactions via the enzyme-substrate binding probability P (j) (k, α).We further assume that shape function is a number, l α and that l k = k.Finally we will drop the dependence of the specificity on α and j, and take it to be a scalar σ.To fix our model, we will take the distortion energy that appears in (3) to be the linear form corresponding to the elastic distortion model [40], do not pose any computational difficulties, and we see that the results of our optimization remain qualitatively unchanged.
These restrictions significantly reduce the dimension of the optimization search, so much so that in certain limits, we can solve the problem analytically 3 .This helps us obtain some useful heuristics (Appendix E) on how to tune the parameters, e.g.N E , N C , σ, and others, in order to generate glycan distributions c of a given complexity.These heuristics inform our more detailed optimization using "realistic" target distributions.
The calculations in Appendix E imply, as one might expect, that the synthesis model needs to be more elaborate, i.e., needs a larger number of cisternae N C or a larger number of enzymes N E , in order to produce a more complex glycan distribution.For a real cell type in a niche, the specific elaboration of the synthesis machinery, would depend on a variety of control costs associated with increasing N E and N C .While an increase in the number of enzymes would involve genetic and transcriptional costs, the costs involved in increasing the number of cisternae could be rather subtle.
Notwithstanding the relative control costs of increasing N E and N C , it is clear from the special case, that increasing the number of cisternae achieves the goal of obtaining an accurate representation of the target distribution.Let us assume that the target distribution is c * k = δ(k − M ) for a fixed M 1, i.e. c * k = 1 when k = M , and 0 otherwise, and that the N E enzymes that catalyse the reactions are highly specific.In this limit, Optimization A reduces to a simple enumeration exercise [46]: clearly one needs N E = M , with one enzyme species for each of the k = 1, . . ., M reactions, in order to generate Pc M .For a single Golgi cisterna with a finite cisternal residence time (finite µ), the chemical synthesis network will generate a significant steady state concentration of lower index glycans Pc k with k < M , contributing to a low fidelity.To obtain high fidelity, one needs multiple Golgi cisternae with a specific enzyme partitioning (E 1 , E 2 , . . ., E M ) with E j enzymes in cisterna j = 1, . . ., N c .This argument can be generalised to the case where the target distribution is a finite sum of delta-functions.The more general case, where the enzymes are allowed to have variable specificity, needs a more detailed study, to which we turn to below.

A. Target distribution from coarse-grained MSMS
As discussed in Sect.II, we obtain the target glycan distribution from glycan profiles for real cells obtained using Mass Spectrometry coupled with determination of molecular structure (MSMS) measurements [22].The raw MSMS data, however, is not suitable as a target distribution.This is because it is very noisy, with chemical noise in the sample and Poisson noise associated with detecting discrete events being the most relevant [47].This means that many of the small peaks in the raw data are not part of the signal, and one has to "smoothen" the distribution to remove the impact of noise.
We use MSMS data from human T-cells [22] for our analysis.As discussed in Sect.II, the Gaussian mixture models (GMM) are often used to approximate distributions with a mixed number of modes or peaks [17], or in our setting, a given fixed complexity.Here, we use a variation of the Gaussian mixture models (see Appendix G for details) to create a hierarchy of increasingly complex distributions to approximate the MSMS raw data.Thus, the 3-GMM and 20-GMM approximations represent the low and high complexity benchmarks, respectively.In Appendix G, we show that the likelihood for the glycan distribution of the human T-cell saturates at 20 peaks.Thus, statistically speaking, the human T-cell glycan distribution is accurately approximated by 20 peaks.
This hierarchy allows us to study the trade-off between the complexity of the target distribution and the complexity of the synthesis model needed to generate the distribution as follows.Let T (i) denote the i th -GMM approximation for the human T-cell MSMS data.We sample this target distribution at indices k = 1, . . ., N s , that represent the glycan indices, and then renormalize to obtain the discrete distribution {T k denote the entropy [16] of the i th -GMM approximation.H(T (i) ) quantifies statistical information in the target distribution T (i) .We evaluate the fidelity of the distribution generated by the synthesis model to this target distribution by the ratio of the Kullback-Leibler distance to the entropy of the target distribution: This normalization allows us evaluate the fidelity of the synthesis model to the target distribution T (i) as a fraction of the total statistical information in the target distribution T (i) .

B. Tradeoffs between number of enzymes, number of cisternae and enzyme specificity to achieve given complexity
We are now in a position to catalogue the main results that follow from an optimization of the parameters of the glycan synthesis machinery to a given target distribution, Figs.provided the sensitivity σ is tightly controlled at σ min (N C , N E , c * ), and there is larger value of (N E , N C ) such that the target can be synthesized even if the control on σ is less tight.
Ungar et al. [26] optimize incoming glycan ratio, transport rate and effective reaction rates in order to synthesize a narrow target distribution centred around a desired glycan.The ability to produce specific glycans without much heterogeneity is an important goal in pharma industry.They define heterogeneity as the total number of glycans synthesized, and show that increasing the number of compartments N C decreases heterogeneity, and increases the concentration of the specific glycan.They also show that changing transport rate does not affect the heterogeneity.
Our results are entirely consistent with theirs -we have shown that D decreases as we increase N C .Thus, if the target distribution has a single sharp peak, increasing N C will reduce the heterogeneity in the distribution.

C. Optimal partitioning of enzymes in cisternae
Having studied the optimum N E , N C , σ to attain a given target distribution with high fidelity, we ask what is the optimal partitioning of the N E enzymes in these N C cisternae?Answering this within our chemical reaction model (Sect.IV A) requires some care, since it incorporates the following enzymatic features: (a) enzymes with a finite specificity σ can catalyse several reactions, although with an efficiency that varies with both the substrate index k and cisternal index j, and (b) every enzyme appears in each cisternae; however their reaction efficiencies depend on  the enzyme levels, the enzymatic reaction rates and the enzyme matching function L, all of which depend on the cisternal index j.
Thus, rather than determining the cisternal partitioning of enzymes, we instead identify chemical reactions that occur with high propensity in each cisternae.For this we define an effective reaction rate R(j, k) for Pc k → Pc k+1 in the .This has implications for glycoengineering (see text) where the task is to produce a particular glycan profile with low heterogeneity [26,46].
According to our model presented in Sect.IV A, the list of reactions with high effective reaction rates in each cisterna, corresponds to a cisternal partitioning of the perfect enzymes.In a future study, we will consider a Boolean version of a more complex chemical model, to address more clearly, the optimal enzyme partitioning amongst cisternae.
Figure 5 (a) (i) shows the heat map of the effective reaction rates in each cisterna for the optimal N E , N C , σ that minimises the normalised KL-distance to the 20 GMM target distribution T (20) (see Fig. 5 (a) (ii)).The optimized glycan profile displayed in Fig. 5 (a) (iii) is very close to the target.An interesting observation from Fig. 5a(i) is that the same reaction can occur in multiple cisternae.
Keeping everything else fixed at the optimal value, we ask whether simply repartitioning the optimal enzymes amongst the cisternae, alters the displayed glycan distribution.In Fig. 5 (b) (i), we have exchanged the enzymes of the fourth and second cisterna.The glycan profile after enzyme partitioning (see Fig. .Thus one may achieve a different glycan distribution by repartitioning enzymes amongst the same number of cisternae [46].

VII. STRATEGIES TO ACHIEVE HIGH GLYCAN DIVERSITY
So far we have studied how the complexity of the target glycan distribution places constraints on the evolution of Golgi cisternal number and enzyme specificity.We now take up another issue, namely, how the physical properties of the Golgi cisternae, namely cisternal number and inter-cisternal transport rate, may drive diversification of glycans [48,49].There is substantial correlative evidence to support the idea that cell types that carry out extensive glycan processing employ larger numbers of Golgi cisternae.For example, the salivary Brunners gland cells secrete mucous that contains heavily O-glycosylated mucin as its major component [50].The Golgi complex in these specialized cells contain 9 − 11 cisternae per stack.Additionally, several organisms such as plants and algae secrete a rather diverse repertoire of large, complex glycosylated proteins, for a variety of functions [51][52][53][54][55][56][57][58][59][60].These organisms possess enlarged Golgi complexes with multiple cisternae per stack [61][62][63][64][65].
In this section, we study how changing the physical parameters in our chemical synthesis model can lead to changes in the diversity of glycan distributions.
We define diversity as the total number of glycan species produced above a specified threshold abundance c th .This last condition is necessary because very small peaks will not be distinguishable in the presence of noise.In computing the diversity from our chemical synthesis model, we have chosen the threshold to be c th = 1/N s , where N s is the total number of glycan species.We have checked that the qualitative results do not depend on this choice, Fig. A6.
Using the sigmoid function (1 + e −x/τ ) −1 as a continuous approximation to the Heaviside function Θ(x), we define the following optimization to achieve the maximal diversity for a given set of parameter values, N E , N C , σ, where, as before, (µ max , µ min ) = (1, 0.01)/min, and (R max , R min ) = (20, 0.018)/min, and c th = 1/N s is the threshold.See Appendix F for details on the parameter estimation.
The results displayed in Fig. 6(a), show that for a fixed specificity σ, the diversity at first increases with the number of cisternae N C , and then saturates at a value that depends on σ.For very high specificity enzymes, one can achieve very high diversity by appropriately increasing N C .This establishes the link between glycan diversity and cisternal number.However, this link is correlative at best, since there are many ways to achieve high glycan diversity -notably by increasing the number of enzymes.
On the other hand, one of the goals of glycoengineering is to produce a particular glycan profile with low heterogeneity [26,46].For low specificity enzymes, the diversity remains unchanged upon increasing the cisternal residence time.For enzymes with high specificity, the diversity typically shows a non-monotonic variation with the cisternal residence time.At small cisternal residence time, the diversity decreases from the peak because of early exit of incomplete oligomers.At large cisternal residence time the diversity again decreases as more reactions are taken to completion.Note that the peak is generally very flat, this is consistent with the results of [26].To get a sharper peak, as advocated for instance by [46], one might need to increase the number of high specificity enzymes N E further.

VIII. DISCUSSION
The precision of the stereochemistry and enzymatic kinetics of these N-glycosylation reactions [2], has inspired a number of mathematical models [23][24][25] that predict the N-glycan distribution based on the activities and levels of processing enzymes distributed in the Golgi-cisternae of mammalian cells, and compare these predictions with N-glycan mass spectrum data.Models such as the KB2005 model [23][24][25] are extremely elaborate (with a network of 22, 871 chemical reactions and 7565 oligosaccharide structures) and require many chemical input parameters.These models have an important practical role to play, that of being able to predict the impact of the various chemical parameters on the glycan distribution, and to evaluate appropriate metabolic strategies to recover the original glycoprofile.Additionally, a recent study by Ungar and coworkers [26,27] shows how physical parameters, such as overall Golgi transit time and cisternal number, can be tuned to engineer a homogeneous glycan distribution.Overall, such models may help predict glycosylation patterns and direct glycoengineering projects to optimize glycoform distributions.
In this paper, we have been interested in the role of glycans as a marker or molecular code of cell identity [4,7,11].
In particular, we have studied one aspect of molecular coding, namely the fidelity of this glycan code generated by enzymatic and transport processes located in the secretory apparatus of the cell.This involves a method of analysis that draws on many different fields, and so it might be useful to provide a short summary of the assumptions, methods and results of the paper: 1.The distribution of glycans at the cell surface is a marker of cell-type identity [2,4,7,11].This glycan distribution can be very complex; it is believed that there is an evolutionary drive for having glycan distributions of high complexity arising from the following considerations, (a) Reliable cell type identification amongst a large set of different cell types in a complex organism, the preservation and diversification of "self-recognition" [5].(b) Consequence of pathogen-mediated selection pressures [2,4,6].(c) Consequence of herd immunity within a heterogenous population of cells of a community [15] or within a single organism [5].
2. The glycans at the cell surface are the end product of a sequential chemical processing via a set of enzymes resident in the Golgi cisternae, and transport across cisternae [4,10,11].Using a fairly general and tractable model for chemical synthesis and transport, we compute the synthesized glycan distribution at the cell surface.Parameters of our synthesis model include the number of enzymes N E , specificity of enzymes σ, number of cisternae N C and transport rates µ.
3. We measure the reliability or fidelity of cell identity [10,11,66] in terms of the error between synthesized glycan distribution on the cell surface from the its internal "target" distribution using the Kullback-Leibler distance D KL [16,17].In our numerical study, we obtain the target distribution for the given cell type by suitable coarse-graining of the MSMS data for the human T-cells [22].We solve a constrained optimization problem for minimising D KL , and study the tradeoffs between N E , N C and σ.
4. The results of the optimization to achieve a given target complexity are summarised in Figs.3-4.Here, we highlight some its direct consequences: (a) Keeping the number of enzymes fixed, a more elaborate transport mechanism (via control of N C and µ) is essential for synthesising high complexity target distributions (Figs.4a, 4b).Fewer cisternae cannot be compensated for by optimising the enzymatic synthesis (via control of parameters R, L and σ).(b) Thus, our study suggests that fidelity of the glycan code generation provides a functional control of Golgi cisternal number.It also provides a quantitative argument for the evolutionary requirement of multiplecompartments, by demonstrating that the fidelity and sensitivity of the glycan code arising from a chemical synthesis that involves multiple cisternae is higher than one that involves a single cisterna (keeping everything else fixed) (Figs.4a, 4b, 4e, 4f).This feature that with multiple cisternae and precise enzyme partitioning, one may generically achieve a highly accurate representation of the target distribution, has been highlighted in an algorithmic model of glycan synthesis [46].(c) Our study shows that for a fixed N C and N E , there is an optimal enzyme specificity that achieves the lowest distance from a given target distribution.As we see in Fig. 4d, this optimal enzyme specificity can be very high for highly complex target distributions.(d) Organisms such as plants and algae, have a diverse repertoire of glycans that are utilised in a variety of functions [51][52][53][54][55][56][57][58][59][60].Our study shows that it is optimal to use low specificity enzymes to synthesize target distributions with high diversity (Fig. 6).However, this compromises on the complexity of the glycan distribution, revealing a tension between complexity and diversity.One way of relieving this tension is to have larger N E and N C .(e) Consider a situation where the environment, and hence the target glycan distribution, fluctuates rapidly.
When synthesis parameters cannot change rapidly enough to track the environment, high specificity enzymes can lead to a lowering of the cell's fitness [67,68].Having slightly sloppy enzymes may give the best selective advantage in a time varying environment.This compromise, between robustness in a changing environment and the demand for complexity, is achieved by having sloppy enzymes, that allows the system to be more evolvable [67,68].However, sloppy enzymes are subject to errors from synthesising the wrong reaction products.In this case, error correcting mechanisms must be in place to ensure fidelity of the glycan code.We leave the role of intra-cellular transport in providing non-equilibrium proof-reading mechanisms to reduce such coding errors, and its optimal adaptive strategies and plasticity in a time varying environment, as a task for the future.
Admittedly the chemical network that we have considered here is much simpler than the chemical network associated with all possible protein modifications in the secretory pathway.For instance, typical N-glycosylation pathways would involve the glycosylation of a variety of GBPs.Further, apart from N-glycosylation, there are other glycoprotein, proteoglycan and glycolipid synthesis pathways [1,2,11].We believe our analysis is generalisable and that the qualitative results we have arrived at would still hold.
To conclude, our work establishes the link between the cisternal machinery (chemical and transport) and optimal coding.We find that the pressure to achieve the target glycan code for a given cell type, places strong constraints on the cisternal number and enzyme specificity [18].An important implication is that a description of the nonequilibrium self-assembly of a fixed number of Golgi cisternae must combine the dynamics of chemical processing and membrane dynamics involving fission, fusion and transport [18][19][20].We believe this is a promising direction for future research.is equivalent to (8).Since v is explicitly known as a function of (R, L), optimization B (B4) is a more tractable optimization problem than (8).Note that in this setting, the function D(σ, N E , N C , c * ) (B4) is independent of the rates µ.
While this optimization is easy to implement, we note that the parameters (e.g., reaction rates, specificity) are not constrained to take only physically relevant values; a legitimate concern is that the absence of such physico-chemical constraints might drive this optimization to physically unrealistic solutions.
There are two possible ways to impose these parameter constraints.One is to impose constraints on the "microscopic" chemical parameters, such as the rate of individual reactions R(j, k, α) and the inter-cisternal transport rate µ (j) .These take into consideration constraints arising from molecular enzymatic processes.The other is to impose constraints on "global" physical parameters, such as the total transport time across the Golgi cisternae and the average enzymatic reaction time.Here, we impose constraints on the microscopic reaction and transport parameters.
The upper and lower bounds on the rates R and µ are estimated in Appendix F : µ max = 1/min (resp.µ min = .01/min)and R max = 20/min (resp.R min = .018/min).and dc for j = 2, . . ., N C .In (D1) and (D2), where the indicator function I(•) is equal to 1 if the argument is true, and zero otherwise and R (j, k, α) is the derivative of R(j, k, α) with respect to k.

Define a vector function
k , . . .c ]. Then (D1) and (D2) can be written as: where the matrix M (k) is given by with The functions A (j) (k) and B (j) (k) involve absolute value and indicator functions; therefore, the differential equation has to be solved in a piecewise manner assuming continuity of solution C(k).
The general solution of (D4) is written in terms of the Magnus Function Ω(k) = ∞ n=1 Ω(n, k), obtained from the Baker-Campbell-Hausdorff formula [69], is the commutator, and the higher order terms in . . .contain higher order nested commutators.
Here, we establish conditions under which the series ∞ n=1 Ω(n, k) that defines solution C(k) to the differential equation (D4) converges.We also solve (D4) for some special cases.
The commutator where The general form of Ω(n, k) is given by [69] Ω where (p We can bound all the matrix elements of Ω(n, k) in the following way The matrix where a lm = S lm (n) Ān for appropriately defined polynomials S l,m (n).Thus, it follows that Since the parameters µ, σ, R(j, k, α), l (j) α and N E are finite and positive, and R (j, k, α) is finite, Ā has a finite upper bound, implying k is always greater than zero, and the series has finite radius of convergence.
While in principle we can obtain the glycan profile for any N E and N C with arbitrary accuracy, assuming R(j, k, α) = R (j) α , we provide explicit formulae for a few representative cases : A representative concentration profile is plotted in Fig. A2(a).The concentration profile consists of two distinct components: an initial exponential decay, and then an exponential rise and fall concentrated around l.The relative weight of these two components is controlled by the sensitivity σ and the rate R.Such explicit formulae can be obtained for any N E > 1, as long as N C = 1.
Appendix E: Capability of the chemical network model to generate complex distributions

Is our glycan synthesis model capable of generating concentration distributions of arbitrary complexity?
In what way do we need to change the parameters N E , N C , σ, . .., in order to generate glycan distributions c of a given complexity?
The purpose of this section, is to obtain some heuristics for this task.
We show in Appendix D that (B4) can be solved analytically in the limit N s 1, because in this limit the glycan index k can be approximated by a continuous variable, and the recursion relations for the steady state glycan concentrations ( 5)-( 6) can be cast as a matrix differential equation.This allows us to obtain an explicit expression for the steady state concentration in terms of the parameters (R, L).
We derive our heuristics from a semi-analytical treatment in the limit N s 1 (Appendix D), which apart from being simple to implement in general, provides an explicit formula for c k for the case N E = N C = 1 (D10).Figures A3(a)-(d) show the glycan profile c k vs. k as one varies the enzyme specificity σ, the reaction rates R and transport rates µ, for two different values of N E and N C .The results in the plots lead us to the following general observations: • Very low specificity enzymes cannot generate complex glycan distributions.Keeping everything else fixed, intermediate or high specificity enzymes can generate glycan distributions of higher complexity by increasing N E or N C (Figs. A3(a),(c)).
• Decreasing the specificity σ or increasing the rates R increases the proportion of higher index glycans.Keeping everything else fixed, changes in the rate R have a stronger impact on the relative weights of the higher index glycans to lower index glycans.The relative weight of the higher index glycans increases with increasing N E and N C (Figs. A3(b)-(d)).
• Keeping everything else fixed, decreasing enzyme specificity increases the spread of the distribution around the peaks (Figs.A3(a),(c)).

Appendix F: Parameter estimation
The typical transport time of glycoproteins across the Golgi complex is estimated to be in the range 15-20 mins.[23], which corresponds to the transport rate, µ = .18/min.We bound the transport rate for our optimization between 0.01/min and 1/min.
Next, we estimate the range of values for the chemical reaction rates.The injection rate q is in the range 100 − 1500 pmol/10 6 cell 24 h [23,24].For our calculation we set q = 387.30pmol/10 6 cells 24 hr = 0.27 pmol/10 6 cells min, the geometric mean of 100 and 1500.We set the range for the enzymatic rate R to be where K (α) M and V (α) max denote the Michaelis constants and V max of the α-th enzyme.The conversion from 1 pmoles/10 6 cells to concentration can be obtained by taking cisternal volume (ν) to be 2.5µm 3 [23,24].This gives 1 pmoles/10 6 cells = 10 −12 moles 10 6 × 2.5 × 10 −18 × 10 3 litre = 400µM.(F1) In Table I we report the parameters for the 8 enzymes taken from Table 3 in [23].From these parameters it follows that 2 , l 1 , l 2 ] = [10,30,50,70]).As σ increases, the distribution becomes more complex -from a single peaked distribution at low σ to a maximum of four-peaked distribution at high σ.The peaks gets sharper, and more well defined as σ increases.Next, we interpret the counts as samples from a "spatial" distribution f .We approximate this distribution as a Gaussian mixture, i.e. f (x) = N =1 γ η(x | µ , σ ), where η(x | µ, σ) denotes the density of a normally distributed random variable with mean µ and standard deviation σ at the location x.In this setting, we assume that each count is a sample from the distribution η(x | µ , σ ) with probability γ .Thus, each count is classified as coming from one of the Gaussian components.adding σ to the optimization vector and then performing the optimization.The sensitivity results (Figs. 4e and 4e) were obtained by approximating the D vs σ graph around σ min with a parabola, the coefficient of the quadratic term being the curvature of the graph at σ min .
A similar numerical scheme was used to optimize diversity.

FIG. 1 .
FIG. 1. Real cells display a complex glycan distribution.(a) Here we take the MSMS data from mouse neutrophils, human neutrophils and human T-cells and approximate these using Gaussian Mixture Models (GMM) of less complexity 3-GMM (left) and more complexity 20-GMM (right).(b) The change in log likelihood with increase in the number of GMM components for mouse and human neutrophils and human T-cells, shows a saturation at large enough values of k, indicating that these glycan distribution are complex.Details appear in Appendix G.

FIG. 2 .
FIG. 2. Enzymatic reaction and transport network in the secretory pathway.Represented here is the array of Golgi cisternae (blue) indexed by j = 1, . . ., NC situated between the ER and PM.Glycan-binding proteins Pc(1) 1

( 1 ) 2 ,
following an MM-reaction (1) catalysed by the appropriate enzyme.The acceptor Pc (1) 2 has the potential of being transformed into Pc

3 - 4 1. 2 .
The normalized KL-distance D(σ, N E , N C , c * ) is a convex function of σ for fixed values for other parameters (Fig.3), i.e. it first decreases with σ and then increases beyond a critical value of σ min .D(σ, N E , N C , c * ) is decreasing in N C and N E for fixed values of the other parameters, and increasing in the complexity of c * for fixed (σ, N C ).The marginal contribution of N C and N E in reducing the normalised distance D is approximately equal (Figs.4a, 4b).The lower complexity distributions can be synthesized with high fidelity with small (N E , N C ), whereas higher complexity distributions require significantly larger (N E , N C ), Figs.4a, 4b.For a typical mammalian cell, the number of enzymes in the N-glycosylation pathway are in the range N E = 10−20[23- 25, 27], Fig.4bwould then suggest that the optimal cisternal number would range from N C = 3 − 8[18].The optimal enzyme specificity σ min (c * , N C ) = argmin σ { D(σ, NE , N C , c * )}, that minimises the error as function of (N C , c * ) with N E fixed at NE , is an increasing function of N C and the complexity of the target distribution c * (Figs.3a, 3b, 4c, 4d).This is consistent with the results in Appendix E where we established that the width of the synthesized distribution is inversely dependent on the specificity σ: since a GMM approximation with fewer peaks has wider peaks, σ min is low, and vice versa.Similar results hold when N C is fixed at NC , and N E is varied (Figs.3c, 3d, 4c, 4d). 3. Let σ min (N C , N E , c * ) denote the value of σ that minimizes D(σ, N C , N E , c * ).Then the second-derivative ∇ 2 σmin D(N C , N E , c * ) = d 2 dσ 2 D(σ, N C , N E , c * ) | σ=σmin denotes the curvature at σ min , and is measure of the sensitivity of D(σ, N E , N C , c * ) to σ for values close to σ min (N E , N C , c * ).∇ 2 σmin D(N C , N E , c * ) is a decreasing function of N C (resp.N E ) for fixed values of (N E , c * ) (resp.(N C , c * )), see Figs. 3, 4e, 4f.Thus, for any target distribution c * there is a minimal value of (N E , N C ) such that the target can be synthesized with high fidelity FIG. 3. Tradeoffs amongst the glycan synthesis parameters, enzyme specificity σ, cisternal number NC and enzyme number NE, to achieve a complex target distribution c * ).(a)-(b) Normalised Kullback-Leibler distance D(σ, NE, NC , c * ) as function of σ and NC (for fixed NE = 3), (c)-(d) D(σ, NE, NC , c * ) as function of σ and NE (for fixed NC = 3), with the target distribution c * set to the 3-GMM (less complex) and 20-GMM (more complex) approximations for the human T-cell MSMS data.D(σ, NE, NC , c * ) is a convex function of σ for each (NE, NC , c * ), decreasing in NC , NE for each (σ, c * ), increasing in the complexity of c * for fixed (σ, NE, NC ).The specificity σmin(c * , NE, NC ) = argmin σ { D(σ, NE, NC , c * )} that minimises the error for given (NE, NC , c * ) is an increasing function of NC , NE and the complexity of the target distribution c * .Furthermore, the curvature of D(σ, NE, NC , c * ) at σmin(NE, NC , c * ), related to sensitivity, is a decreasing function of NC , NE.

( a )FIG. 4 .
FIG. 4. Fidelity of glycan distribution and optimal enzyme properties to achieve a complex target distribution.The target c * is taken from 3-GMM (less complex) and 20-GMM (more complex) approximations of the human T-cell MSMS data.(a)-(b) Minimum normalised KL divergence minσ{ D(σ, NC , NE, c * )} as a function of (NE, NC ).More complex distributions require either a larger value NE or NC .The marginal impact of increasing NE and NC on D is approximately equal.(c)-(d) Optimum enzyme specificity σmin obtained from minσ{ D(σ, NC , NE, c * )} as a function of (NE, NC ).σmin increases with increasing NE or NC .To synthesize the more complex 20 GMM approximation with high fidelity requires enzymes with higher specificity σmin compared to those needed to synthesize the broader, less complex 3-GMM approximation.(e) -(f) Sensitivity ln d 2 dσ 2 D σ min of the normalised Kullback-Leibler distance D(σ, NC , NE, c * ) as a function of (NE, NC ).Increasing NE or NC decreases this sensitivity implying the specificity σ does not need to be tuned very carefully if NE, NC are high.

FIG. 5 .
FIG. 5. Optimal enzyme partitioning in cisternae.(a) Heat map of the effective reaction rates in each cisterna (representing the optimal enzyme partitioning) and the steady state concentration in the last compartment (c (N C ) ) for the 20-GMM target distribution.Here NE = 5, NC = 7, normalised DKL(T (20) c (N C ) )/H(T (20) ) = 0.11.(b) Effective Reaction rates after swapping the optimal enzymes of the fourth and second cisternae.The displayed glycan profile is considerably altered from the original profile.

FIG. 6 .
FIG. 6. Strategies for achieving high glycan diversity.Diversity versus NC and transport rate µ at various values of specificity σ for fixed NE = 3.(a) Diversity vs. NC at optimal transport rate µ.Diversity initially increases with NC , but eventually levels off.The levelling off starts at a higher NC when σ is increased.These curves are bounded by the σ = 0 curve.(b) Diversity vs. cisternal residence time (µ −1 ) in units of the reaction time (R −1 min ) at various value of σ, for fixed NC = 4 and NE = 10.This has implications for glycoengineering (see text) where the task is to produce a particular glycan profile with low heterogeneity[26,46].
5 (b) ((iii)) is now completely altered (compare Fig. 5 (b) (ii) with Fig. 5 (b) (iii)) FIG. A1.Flow chart showing the optimization schemes for Optimization A and B. We prove that D A min = D B min by showing the set of all c is equal to the set of all v.We additionally establish that the optimum vmin = cmin.

= 5 ×
FIG.A4.The binned MS data (blue) approximates the raw MS data (red) very well.We use this binned data for GMM approximation of the MS data.
FIG. A5.Log likelihood vs. number of components (N ) in the GMM.We see that the log likelihood saturates at around N = 20, thus 20-GMM is a very good representation of the MS-data from human T-cells.The different symbols are for (a) different values of the maximum intensity Imax = 50, 100, 200 and (b) different values of the number of i.i.d samples Ni = 500, 1000, 2000, showing the insensitivity of the log likelihood to the value of Imax and Ni.

(a) c th = 1
FIG. A6.Diversity vs. NC for different values of σ keeping NE = 1 fixed, for three different values of the threshold, c th = 1 Ns , 1 2Ns , 1 4Ns .Changing the value of the threshold c th , only changes the saturation value of the diversity curve.