Glycan processing in the Golgi as optimal information coding that constrains cisternal number and enzyme specificity

Version of Record

Accepted for publication after peer review and revision.

Download
Cite
Share
CommentOpen annotations (there are currently 0 annotations on this page).

Version of Record published: May 31, 2022 (This version)
Accepted Manuscript published: February 17, 2022 (Go to version)
Accepted: January 31, 2022
Received: January 4, 2022
Preprint posted: May 19, 2020 (Go to version)

1. Of interest
Substrate stiffness impacts early biofilm formation by modulating Pseudomonas aeruginosa twitching motility

Sofia Gomez, Lionel Bureau ... Sigolene Lecuyer

Research Article Updated Jun 9, 2023
Further reading

Article
Figures and data
Abstract
Editor's evaluation
Introduction
Model
Results
Discussion
Appendix 1
Appendix 2
Appendix 3
Appendix 4
Appendix 5
Appendix 6
Appendix 7
Appendix 8
Appendix 9
Appendix 10
Appendix 11
Data availability
References
Decision letter
Author response
Article and author information
Metrics

Abstract

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 realizes the desired target glycan distribution for a particular cell type and niche. In this article, we construct a simplified chemical synthesis model to quantitatively analyse the trade-offs between the number of cisternae, and the number and specificity of enzymes, required to synthesize a prescribed target glycan distribution of a certain complexity to within a given fidelity. 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 a fixed number of enzymes and cisternae, there is an optimal level of specificity (promiscuity) of enzymes that achieves the target distribution with high fidelity. The geometry of the fidelity landscape in the multidimensional space of the number and specificity of enzymes, inter-cisternal transfer rates, and number of cisternae provides a measure for robustness and identifies stiff and sloppy directions. Our results show how the complexity of the target glycan distribution and number of glycosylation enzymes places functional constraints on the Golgi cisternal number and enzyme specificity.

Editor's evaluation

This article contributes to an important and largely unexplored topic in cell biology: the understanding of glycosylation. The authors introduce a mathematical model of glycosylation in the Golgi apparatus and use the model to investigate how the complexity (diversity) and fidelity of the plasma membrane glycan distribution depend on parameters such as the number of Golgi cisternae or enzyme specificity. The article is well written and makes the effort to present a rather complex topic in an accessible way by leaving some of the details in the appendices.

https://doi.org/10.7554/eLife.76757.sa0

Introduction

A majority of the proteins synthesized in the endoplasmic reticulum (ER) are transferred to the Golgi cisternae for further chemical modification by glycosylation (Alberts, 2002), 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 (Varki, 2009; Cummings and Pierce, 2014). 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 (Varki, 2009; Cummings and Pierce, 2014; Varki, 2017; Drickamer and Taylor, 1998; Gagneux and Varki, 1999). This glycan code (Gabius, 2018; Dwek, 1996), representing information (Winterburn and Phelps, 1972) about the cell, is generated dynamically, following the biochemistry of sequential enzymatic reactions and the biophysics of secretory transport (Varki, 2017; Varki, 1998; Pothukuchi et al., 2019).

In this article, 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 (Gabius, 2018; Varki, 2017; Pothukuchi et al., 2019). While the functional consequences of glycan alterations have been well studied, the glycan code has remained an enigma (Gabius, 2018; Pothukuchi et al., 2019; Bard and Chia, 2016; D’Angelo et al., 2013). We study the fidelity of molecular code generation, that is, the precision and reliability with which the glycan distribution is created. While it has been recognized that fidelity of the glycan code is necessary for reliable cellular recognition (Demetriou et al., 2001), a quantitative measure of fidelity of the mechanism and the constraints that fidelity requirements put on cellular structure and organization are lacking.

There are two aspects of the cell-type-specific glycan code and the code generation mechanism that have an important bearing on quantifying fidelity. The first is that extant glycan distributions have high complexity (section ‘Complexity of glycan code’), 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’ (Drickamer and Taylor, 1998), (b) pathogen-mediated selection pressures (Varki, 2009; Varki, 2017; Gagneux and Varki, 1999), and (c) herd immunity within a heterogenous population of cells of a community (Wills and Green, 1995) or within a single organism (Drickamer and Taylor, 1998). We interpret this to mean that the target distribution of glycans for a given cell type is complex; in section ‘Complexity of glycan code’, we define a quantitative measure for complexity and demonstrate its implications. 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 (Varki, 2017; Varki, 1998; Pothukuchi et al., 2019); 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 section ‘Synthesis of glycans in the Golgi cisternae’.

Here, we define fidelity as the minimum achievable Kullback-Leibler (KL) divergence (Cover and Thomas, 2012; MacKay, 2003) between the synthesized distribution of glycans and the target glycan distribution as a function of given cellular parameters, such as the number and specificity of enzymes, inter-cisternal transfer rates, and number of cisternae (section ‘Optimization problem’). Using a simplified chemical synthesis model, we analyse the trade-offs between the number of cisternae and the number and specificity of enzymes in order to achieve a prescribed target glycan distribution with high fidelity (section ‘Results of optimization’). Our analysis leads to a number of interesting results, a few of which we list here:

First, since an important function of the glycan spectrum is cell type/niche identification, it seems natural to relate glycan complexity to organismal complexity taken to be associated with the number of cell types in the organism (Carroll, 2001; Bonner, 1998). Here, we provide a measure of the complexity of the glycan distribution of a given cell type using mass spectrometry coupled with determination of molecular structure (MSMS) data. Using this we have analysed the MSMS data from hydra, planaria, and mammalian cells. We find that the complexity of the glycan distribution indeed correlates with the organism complexity.
Constructing a high-fidelity representation of a complex target distribution, such as those observed in real cells, requires a complex Golgi machinery with multiple cisternae, precise enzyme partitioning, and control on enzyme specificity. This definition of fidelity of the glycan code allows us to provide a quantitative argument for the evolutionary requirement of multiple compartments. While it is possible to produce complex glycan distributions in one compartment using a large number of enzymes, such a design would inevitably require a more elaborate genetic cost.
Within our synthesis model, an increase in the number of Golgi cisternae drives an increase in the glycan complexity, keeping everything else fixed.
We explore the geometry of the fidelity landscape in the multidimensional space of the number and specificity of enzymes, inter-cisternal transfer rates, and number of cisternae. This allows us to discuss issues such as robustness to noise, and stiff and sloppy directions in this multidimensional space.
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. 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. Promiscuous enzymes bring in the potential for evolvability (Kirschner and Gerhart, 2008); promiscuity allows the system to be stable to random mutations in proteins or variations in the target distribution.

Thus, our results imply that the pressure to produce the target glycan code for a given cell type with high fidelity places strong constraints on the cisternal number and enzyme specificity (Sengupta and Linstedt, 2011). Taken together, our quantitative analysis of the trade-offs has deep implications for the non-equilibrium self-assembly of Golgi cisternae and suggests that the control of cisternal number must involve a coupling of non-equilibrium self-assembly of cisternae with enzymatic chemical reaction kinetics (Glick and Malhotra, 1998). This combined dynamics of chemical processing with non-equilibrium membrane dynamics involving fission, fusion, and transport (Sachdeva et al., 2016; Sens and Rao, 2013) opens up a new direction for future research.

Complexity of glycan code

Since each cell type (in a niche) is identified with a distinct glycan profile (Gabius, 2018; Varki, 2017; Pothukuchi et al., 2019), and this glycan profile is noisy because of the stochastic noise associated with the synthesis and transport (Pothukuchi et al., 2019; Bard and Chia, 2016; D’Angelo et al., 2013), 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. Our task is to identify a quantitative measure for the complexity of a glycan profile such that a set of more complex glycan profiles is able to support a larger number of well-separated profiles, and therefore, a larger number of cell types, or equivalently, a more complex organism (a rigorous definition of complexity can be given in terms of the KL metric [Cover and Thomas, 2012; MacKay, 2003] between two glycan profiles. We declare that two profiles are distinguishable only if the KL distance between the profiles is more than a given tolerance. This tolerance is an increasing function of the noise. We define the complexity of a set of possible glycan profiles as the size of the largest subset such that the KL distance of any pair of profiles is larger than the tolerance). Furthermore, we would like to be able to estimate the complexity of a glycan profile from molecular structure (MSMS) measurements (Cummings and Crocker, 2020; Subramanian et al., 2018; Sahadevan et al., 2014).

In order to identify such a quantitative measure of complexity, we first need a consistent way of smoothening or coarse-graining the raw glycan profiles obtained from MSMS measurements to remove measurement and synthesis noise. Here, we denoise the glycan profile by approximating it by a Gaussian mixture model (GMM) with a specified number of components that are supported on a finite set of indices (Bacharoglou, 2010). Since the size of the set of all possible $m$ -component Gaussian densities is an increasing function of $m$ , we define the complexity of a mixture of Gaussians as the number of components $m$ . Figure 1 demonstrates that the value of $m$ at which the $m$ -component GMM approximation of the target profile saturates is a good measure of complexity. Using this definition we see that the complexity of the glycan profiles of various organisms correlates well with the number of cell types in an organism (details of the procedure are given in Appendix 1). We will now describe a general model of the cellular machinery that is capable of synthesizing glycans of any complexity. We expect that cells need a more elaborate mechanism to produce profiles from a more complex set.

Figure 1

Download asset Open asset

Living cells display a complex glycan distribution.

(a) 3-Gaussian mixture model (GMM) and 20-GMM approximation for the relative abundance of glycans taken from mass spectrometry coupled with determination of molecular structure (MSMS) data of planaria *Schmidtea mediterranea*, *Hydra magnipapillata,* and *human* neutrophils. (b) The change in the Kullback–Leibler (KL) divergence $D (p_{T} ∥ p_{G M M}^{(m)})$ as a function of the number of GMM components $m$ . The KL divergence for planaria saturates at $m = 5$ , for hydra at $m = 11$ , and for *human* cells at $m = 20$ . Thus, the number of components required to approximate the glycan profile correlates well with the complexity of the organism. Details are given in Appendix 1.

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 PM, as depicted in Figure 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 (Varki, 2009). 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, for example, UDP-Gal, is transferred to the acceptor substrate, by a Michaelis–Menten (MM)-type reaction (Varki, 2009)

\begin{array}{ll} Acceptor + glycosyl donor + Enzyme \\ ⇌_{ω_{b}}^{ω_{f}} [Acceptor \cdot glycosyl donor \cdot Enzyme] \overset{ω_{c}}{\to} glycosylated acceptor + nucleotide + Enzyme \end{array}

Figure 2

Download asset Open asset

Enzymatic reaction and transport network in the secretory pathway.

Represented here is the array of Golgi cisternae (blue) indexed by $j = 1, \dots, N_{C}$ situated between the endoplasmic reticulum (ER) and plasma membrane (PM). Glycan-binding proteins $P c_{1}^{(1)}$ are injected from the ER to cisterna-1 at rate $q$ . Superimposed on the Golgi cisternae is the transition network of chemical reactions (column) – inter-cisternal transfer (rows), the latter with rates $μ^{(j)}$ . $P c_{1}^{(1)}$ denotes the acceptor substrate in compartment $j$ and the glycosyl donor $c_{0}$ is chemostated in each cisterna. This results in a distribution (relative abundance) of glycans displayed at the PM (red curve), which is representative of the cell type.

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 catalyzed by the enzymes resident in the second cisterna occurs. This chemical processing and inter-cisternal transfer continue until the last cisterna, thereupon the fully processed glycans are displayed at the PM (Varki, 2009). The network of chemical processing and inter-cisternal transfer forms the basis of 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 allowed values for the reaction rates and chemical abundances, we use the elaborate enzymatic reaction models, such as the KB2005 model (Umaña and Bailey, 1997; Krambeck et al., 2009; Krambeck and Betenbaugh, 2005) (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. For the allowed rates of cisternal transfer, we rely on the recent study by Ungar and coworkers (Fisher et al., 2019; Fisher and Ungar, 2016), whose study shows how the overall Golgi transit time and cisternal number can be tuned to engineer a homogeneous glycan distribution.

Model

Chemical reaction and transport network in cisternae

We consider an array of $N_{C}$ Golgi cisternae, labelled by $j = 1, \dots, N_{C}$ , between the ER and PM (Figure 2). GBPs, denoted as $P c_{1}^{(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 jth cisterna is chemostated (Varki, 2009; Hirschberg et al., 1998; Caffaro and Hirschberg, 2006; Berninsone and Hirschberg, 2000), thus in our model we hold its concentration $c_{0}^{(j)}$ constant in time for the jth cisterna. The acceptor $P c_{1}^{(1)}$ reacts with $c_{0}^{(1)}$ to form the glycosylated acceptor $P c_{2}^{(1)}$ , following an MM reaction (1) catalyzed by the appropriate enzyme. The acceptor $P c_{2}^{(1)}$ has the potential of being transformed into $P c_{3}^{(1)}$ , and so on, provided the requisite enzymes are present in that cisterna. This leads to the sequence of enzymatic reactions $P c_{1}^{(1)} \to P c_{2}^{(1)} \to \dots P c_{k}^{(1)} \to \dots$ , where $k$ enumerates the sequence of glycosylated acceptors using a consistent scheme (such as in Umaña and Bailey, 1997). 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 Figure 2. Although, in this article, 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 graph.

Let $N_{s} - 1$ denote the maximum number of possible glycosylation reactions in each cisterna $j$ , catalyzed by enzymes labelled as $E_{α}^{(j)}$ , with $α = 1, \dots, 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 in Figure 2 is $N_{s} \times N_{C}$ . For the N-glycosylation pathway in a typical mammalian cell, $N_{s} = 2 \times 10^{4}$ , $N_{E}$ = 10–20, and $N_{C}$ = 4–8 (Umaña and Bailey, 1997; Krambeck and Betenbaugh, 2005; Krambeck et al., 2009; Fisher and Ungar, 2016). We account for the fact that the enzymes have specific cisternal localization by setting their concentrations to zero in those cisternae where they are not present.

The action of enzyme $E_{α}^{(j)}$ on the substrate $P c_{k}^{(j)}$ in cisterna $j$ is given by

P c_{k}^{(j)} + E_{α}^{(j)} ⇌_{ω_{b} (j, k, α)}^{ω_{f} (j, k, α) c_{0}^{(j)}} [E_{α}^{(j)} - P c_{k}^{(j)} - c_{0}^{(j)}] \overset{ω_{c} (j, k, α)}{\to} P c_{k + 1}^{(j)} + E_{α}^{(j)}

where $k = 1, \dots N_{s} - 1$ . 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 $α$ , which parametrize the MM reactions (Price and Stevens, 1999). For instance, structural studies on glycosyltransferase-mediated synthesis of glycans (Moremen and Haltiwanger, 2019) would suggest that the forward rate $ω_{f}$ depends on the binding energy of the enzyme $E_{α}^{(j)}$ to acceptor substrate $P c_{k}^{(j)}$ and a physical variable that characterizes the cisternae.

A potential candidate for such a cisternal variable is pH (Kellokumpu, 2019), whose value is maintained homeostatically in each cisterna (Casey et al., 2010); changes in pH can affect the shape of an enzyme (substrate) or its charge properties, and in general the reaction efficiency of an enzyme has a pH optimum (Price and Stevens, 1999). Another possible candidate for a cisternal variable is membrane bilayer thickness (Dmitrieff et al., 2013); indeed, both pH (Llopis et al., 1998) and membrane thickness are known to have a gradient across the Golgi cisternae. We take $ω_{f} (j, k, α) \propto P^{(j)} (k, α)$ , where $P^{(j)} (k, α) \in (0, 1)$ is the binding probability of enzyme $E_{α}^{(j)}$ with substrate $P c_{k}^{(j)}$ , and define the binding probability $P^{(j)} (k, α)$ using a biophysical model, similar in spirit to the Monod-Wyman-Changeux model of enzyme kinetics (Monod et al., 1965; Changeux and Edelstein, 2005) that depends on enzyme-substrate-induced fit.

Let $ℓ_{α}^{(j)}$ and $ℓ_{k}$ denote, respectively, the optimal ‘shape’ for enzyme $E_{α}^{(j)}$ and the substrate $P c_{k}^{(j)}$ . We assume that the mismatch (or distortion) energy between the substrate $k$ and enzyme $E_{α}^{(j)}$ is $∥ ℓ_{k} - ℓ_{α}^{(j)} ∥$ , with a binding probability given by

P^{(j)} (k, α) = \exp (- σ_{α}^{(j)} ∥ ℓ_{k} - ℓ_{α}^{(j)} ∥)

where $∥ \cdot ∥$ is a distance metric defined on the space of $ℓ_{α}^{(j)}$ (e.g. the square of the $ℓ_{2}$ -norm would be related to an elastic distortion model [Savir and Tlusty, 2007]) and the vector $σ \equiv [σ_{α}^{(j)}]$ parametrizes enzyme specificity. This distortion model captures the above idea that the reaction between the flexible enzyme and fixed substrate is facilitated by an induced fit. A large value of $σ_{α}^{(j)}$ indicates a highly specific enzyme, a small value of $σ_{α}^{(j)}$ indicates a promiscuous enzyme. It is recognized that the degree of enzyme specificity or sloppiness is an important determinant of glycan distribution (Varki, 2009; Roseman, 2001; Hossler et al., 2007; Yang et al., 2018).

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 (Umaña and Bailey, 1997; Krambeck et al., 2009; Krambeck and Betenbaugh, 2005). Then the usual MM steady-state conditions for (2), which assumes that the concentration of the intermediate enzyme-substrate complex does not change with time, imply that

[E_{α}^{(j)} - P c_{k}^{(j)} - c_{0}^{(j)}] = \frac{ω_{f} (j, k, α) c_{0}^{(j)}}{ω_{b} (j, k, α) + ω_{c} (j, k, α)} E_{α}^{(j)} c_{k}^{(j)} . {}

where $c_{k}^{(j)}$ is the concentration of the acceptor substrate $P c_{k}^{(j)}$ in compartment $j$ .

Together with the constancy of the total enzyme concentration, ${[E_{α}^{(j)}]}_{t o t} = E_{α}^{(j)} + \sum_{k = 1}^{N_{s}}$ $[E_{α}^{(j)} - P c_{k}^{(j)} - c_{0}^{(j)}]$ , this immediately fixes the kinetics of product formation (not including inter-cisternal transport),

\frac{d c_{k + 1}^{(j)}}{d t} = \sum_{α = 1}^{N_{E}} \frac{V (j, k, α) P^{(j)} (k, α) c_{k}^{(j)}}{M (j, k, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(j)} (k^{'}, α) c_{k^{'}}^{(j)}}{M (j, k^{'}, α)})} k = 1, \dots, N_{S}; j = 1, \dots, N_{C}

where

M (j, k, α) = \frac{ω_{b} (j, k, α) + ω_{c} (j, k, α)}{ω_{f} (j, k, α) c_{0}^{(j)}} P^{(j)} (k, α)

and

V (j, k, α) = ω_{c} (j, k, α) {[E_{α}^{(j)}]}_{t o t} .

This reparametrization of the reaction rates $ω_{f}, ω_{b}, ω_{c}$ in terms of $M, V$ is convenient since it relates to experimentally measurable parameters $V_{m a x}$ and MM constant $K_{M}$ , for each $(j, k, α)$ , which can be easily read out (see Appendix 2). As is the usual case, the maximum velocity $V_{m a x}$ is not an intrinsic property of the enzyme because it is dependent on the enzyme concentration ${[E_{α}^{(j)}]}_{t o t}$ ; while $K_{M} (j, k, α) = M (j, k, α) c_{0}^{(j)} / P^{(j)} (k, α)$ is an intrinsic parameter of the enzyme and the enzyme-substrate interaction. The enzyme catalytic efficiency, the so-called “ $k_{c a t} / K_{M}$ ” $\propto P^{(j)} (k, α)$ , is high for perfect enzymes (Bar-Even et al., 2015) 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 3, we display the complete set of equations that describe the changes in the substrate concentrations $c_{k}^{(j)}$ with time. These kinetic equations automatically obey the conservation law for the protein concentration ( $p$ ). At steady state, these kinetic equations lead to a set of nonlinear recursion equations (15)-(16) that are displayed in Appendix 3, which can be solved numerically to obtain the steady-state glycan concentrations, $c \equiv c_{k}^{(j)}$ , as a function of the independent vectors $M \equiv [M (j, k, α)]$ , $V \equiv [V (j, k, α)]$ , and $L \equiv [P^{(j)} (k, α)]$ , the transport rates $μ \equiv [μ^{(j)}]$ and specificity, $σ \equiv [σ_{α}^{(j)}]$ .

Optimization problem

Let $c^{*}$ denote the ‘target’ concentration distribution, normalized to the distribution so that $\sum_{k = 1}^{N_{s}} c_{k}^{*} = 1$ , for a particular cell type, that is, the goal of the sequential synthesis mechanism described in the section ‘Chemical reaction and transport network in cisternae’ is to approximate $c^{*}$ . Let $\bar{c}$ denote the normalized steady-state glycan concentration distribution displayed on the PM. Then Equation 16 implies that ${\bar{c}}_{k} = μ^{(N_{C})} c_{k}^{(N_{C})}$ , $k = 1, \dots, N_{s}$ . We measure the fidelity $F (c^{*} ‖ \bar{c})$ between the $c^{*}$ and $\bar{c}$ by the ratio of the KL divergence $D (c^{*} ‖ c)$ (Cover and Thomas, 2012; MacKay, 2003) to the entropy $H (c^{*})$

F (c^{*} ‖ \bar{c}) := \frac{D (c^{*} ‖ \bar{c})}{H (c^{*})} = \frac{\sum_{k = 1}^{N_{s}} c_{k}^{*} \ln (\frac{c_{k}^{*}}{{\bar{c}}_{k}}) = \sum_{k = 1}^{N_{s}} c_{k}^{*} \ln (\frac{c_{k}^{*}}{c_{k}^{(N_{C})} μ^{(N_{C})}})}{\sum_{k = 1}^{N_{s}} c_{k}^{*} \ln (1 / c_{k}^{*})}

The reason why we divide the KL divergence by the entropy of the target distribution is to enable comparison of the fidelity of the mechanism across target distributions of different complexity. Note that high fidelity corresponds to low values of $F (c^{*} ‖ \bar{c})$ , vice versa.

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

O p t i m i z a t i o n A : \bar{D} (σ, N_{E}, N_{C}, c^{*}) := min_{μ, M, V, L \geq 0} F (c^{*} ‖ \bar{c}),

where we emphasize that the optimum fidelity $\bar{D} (σ, N_{E}, N_{C}, c^{*})$ is a function of $(σ, N_{E}, N_{C}, c^{*})$ . Note that there is a 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 trade-offs and changes of parameters are related to evolutionary timescales.

Optimization A, though well-defined, is a hard problem since the steady-state concentrations (16) are not explicitly known in terms of the parameters $(μ, M, V, L)$ . In Appendix 4, we formulate an alternative problem Optimization B in which the steady-state concentrations are defined explicitly in terms of new parameters $μ, R$ , and $L$ , and 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 5, we describe the variant of the sequential quadratic programming (SQP) (Boyd and Vandenberghe, 2004), which we use to numerically solve the optimization problem.

Results

The dimension of the optimization search space is extremely large $\approx O (N_{s} \times N_{E} \times N_{C})$ . To make the optimization search more manageable, we make the following simplifying assumptions:

We ignore the $k$ -dependence of the vectors $(M, V)$ , or alternatively of $R$ – see Appendix 4 for details.
The enzyme-substrate-binding probability $P^{(j)} (k, α)$ is still dependent on the substrate $k$ . We assume that the shape function is a scalar (a length), that is, $l_{α}^{(j)} = ℓ_{α}^{(j)}$ . It further simplifies the algebra to assume that the lengths of the substrates are integer multiples of a basic unit (which we take to be 1), that is, $ℓ_{k} = k$ . The norm that appears in (3) is taken to be the absolute value difference $| l_{k} - l_{α}^{(j)} |$ . Other metrics, such as ${| l_{k} - l_{α}^{(j)} |}^{2}$ , corresponding to the elastic distortion model (Savir and Tlusty, 2007), do not pose any computational difficulties, and we see that the results of our optimization remain qualitatively unchanged.
We drop the dependence of the specificity on $α$ and $j$ , and take it to be a scalar $σ$ .

These restrictions significantly reduce the dimension of the optimization search, so much so that in certain limits we can solve the problem analytically (in Appendix 6, we show that Equation 21 can be solved analytically in the limit $N_{s} ≫ 1$ since the glycan index $k$ can be approximated by a continuous variable, and the recursion relations for the steady-state glycan concentrations Equations 15–16 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)$ ). This helps us obtain some useful heuristics (Appendix 6) on how to tune the parameters, for example, $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 6 imply, as one might expect, that the synthesis model needs to be more elaborate, that is, 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. Suppose the target distribution $c_{k}^{*} = δ (k - M)$ for a fixed $M ≫ 1$ , that is, $c_{k}^{*} = 1$ when $k = M$ , and 0 otherwise, and that the $N_{E}$ enzymes that catalyze the reactions are highly specific. In this limit, Optimization A reduces to a simple enumeration exercise (Jaiman and Thattai, 2018): clearly, one needs $N_{E} = M$ enzymes, one for each $k = 1, \dots, M$ reactions, in order to generate $P c_{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 $P c_{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}, \dots, E_{M})$ with $E_{j}$ enzymes in cisterna $j = 1, \dots, N_{c}$ . This argument can be generalized 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 next.

Target distribution from coarse-grained MSMS

As discussed in the section ‘Complexity of glycan code’, we obtain the target glycan distribution from glycan profiles for real cells using MSMS measurements (Cummings and Crocker, 2020). 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 (Du et al., 2008). 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 (Cummings and Crocker, 2020) for our analysis. As discussed in the section ‘Complexity of glycan code’, the GMMs are often used to approximate distributions with a mixed number of modes or peaks (MacKay, 2003), or in our setting, a given fixed complexity. Here, we use a variation of the GMMs (see Appendix 1 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 1, we show that the likelihood for the glycan distribution of the human T-cell saturates at 20 peaks. Thus, statistically 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$ -component GMM approximation for the human T-cell MSMS data. We sample this target distribution at indices $k = 1, \dots, N_{s}$ , that represent the glycan indices, and then renormalize to obtain the discrete distribution ${T_{k}^{(i)}, k = 1, \dots, N_{s}}$ . To highlight the role of target distribution complexity, we focus on the 3-GMM $T^{(3)}$ (low complexity) and 20-GMM approximation $T^{(20)}$ (high complexity) in describing our results.

Trade-offs between number of enzymes, number of cisternae, and enzyme specificity to achieve given complexity

We summarize the main results that follow from an optimization of the parameters of the glycan synthesis machinery to a given target distribution in Figure 3 and Figure 4.

Figure 3

Download asset Open asset

Trade-offs amongst the glycan synthesis parameters, enzyme specificity $σ$ , cisternal number $N_{C}$ , and enzyme number $N_{E}$ to achieve a complex target distribution $c^{*}$ .

(**a, b**) Normalised Kullback–Leibler distance $\bar{D} (σ, N_{E}, N_{C}, c^{*})$ as a function of $σ$ and $N_{C}$ (for fixed $N_{E} = 3$ ), (**c, d**) $\bar{D} (σ, N_{E}, N_{C}, c^{*})$ as a function of $σ$ and $N_{E}$ (for fixed $N_{C} = 3$ ), with the target distribution $c^{*}$ set to the 3-Gaussian mixture model (GMM) (less complex) and 20-GMM (more complex) approximations for the *human* T-cell mass spectrometry coupled with determination of molecular structure (MSMS) data. $\bar{D} (σ, N_{E}, N_{C}, c^{*})$ is a convex function of $σ$ for each $(N_{E}, N_{C}, c^{*})$ , decreasing in $N_{C}, N_{E}$ for each $(σ, c^{*})$ , increasing in the complexity of $c^{*}$ for fixed $(σ, N_{E}, N_{C})$ . The specificity $σ_{\min} (c^{*}, N_{E}, N_{C}) = {argmin}_{σ} {\bar{D} (σ, N_{E}, N_{C}, c^{*})}$ that minimizes the error for given $(N_{E}, N_{C}, c^{*})$ is an increasing function of $N_{C}, N_{E}$ and the complexity of the target distribution $c^{*}$ .

Figure 4

Download asset Open asset

Fidelity of glycan distribution and optimal enzyme properties to achieve a complex target distribution.

The target $c^{*}$ is taken from 3-Gaussian mixture model (GMM) (less complex) and 20-GMM (more complex) approximations of the *human* T-cell mass spectrometry coupled with determination of molecular structure (MSMS) data. (**a, b**) Optimum fidelity $\min_{σ} {\bar{D} (σ, N_{C}, N_{E}, c^{*})}$ as a function of $(N_{E}, N_{C})$ . More complex distributions require either a larger $N_{E}$ or $N_{C}$ . The marginal impact of increasing $N_{E}$ and $N_{C}$ on the fidelity $\bar{D}$ is approximately equal. (**c, d**) Enzyme specificity $σ_{\min}$ that achieves $\min_{σ} {\bar{D} (σ, N_{C}, N_{E}, c^{*})}$ as a function of $(N_{E}, N_{C})$ . $σ_{m i n}$ increases with increasing $N_{E}$ or $N_{C}$ . 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.

The optimal fidelity $\bar{D} (σ, N_{E}, N_{C}, c^{*})$ is a convex function of $σ$ for fixed values for other parameters (see Figure 3), that is, it first decreases with $σ$ and then increases beyond a critical value of $σ_{\min}$ .

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})$ (see Figure 4a and b). For a typical mammalian cell, the number of enzymes in the N-glycosylation pathway is in the range $N_{E} = 10 - 20$ (Umaña and Bailey, 1997; Krambeck and Betenbaugh, 2005; Krambeck et al., 2009; Fisher and Ungar, 2016), Figure 4b would then suggest that the optimal cisternal number would range from $N_{C} = 3 - 8$ (Sengupta and Linstedt, 2011).

The fidelity $\bar{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 improving fidelity $\bar{D}$ is approximately equal (see Figure 4a and b). We discuss the origin of this symmetry later in this section.
The optimal enzyme specificity $σ_{\min} (c^{*}, N_{C}) = {argmin}_{σ} {\bar{D} (σ, {\bar{N}}_{E}, N_{C}, c^{*})}$ , which minimizes the error as function of $(N_{C}, c^{*})$ with $N_{E}$ fixed at ${\bar{N}}_{E}$ , is an increasing function of $N_{C}$ and the complexity of the target distribution $c^{*}$ (Figure 3a and b and Figure 4c and d). This is consistent with the results in Appendix 6 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 ${\bar{N}}_{C}$ , and $N_{E}$ is varied (see Figure 3c and d and Figure 4c and d).

Our results are consistent with those in Fisher et al., 2019. They optimize incoming glycan ratio, transport rate, and effective reaction rates in order to synthesize a narrow target distribution centred around the desired glycan. The ability to produce specific glycans without much heterogeneity is an important goal in the pharmaceutical 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 the effect of compartments in reducing heterogeneity cannot be compensated by changing the transport rate. Our results are entirely consistent with theirs – we have shown that $\bar{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.

We insert an important cautionary note here. It would seem that the results in Figure 4 imply that there is an approximate $N_{E} - N_{C}$ symmetry in the model, that is, increasing either $N_{E}$ or $N_{C}$ affects the fidelity, optimal enzyme specificity, and the sensitivity in approximately the same way. This would be an erroneous inference, and is a consequence of the distortion model we have used for calculating the binding probabilities of substrates with enzymes. The root cause for this apparent symmetry is that we have allowed for all enzymes to catalyze reactions in all cisternae (albeit with different efficiencies). This symmetry is violated by simply restricting the activity of the enzymes to be dependent on the cisternae. A simple realization of this in terms of the distortion model is given in Appendix 7.

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 the context of our chemical reaction model (section ‘Chemical reaction and transport network in cisternae’) requires some care since it incorporates the following enzymatic features: (a) enzymes with a finite specificity $σ$ can catalyze 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$ .

Therefore, instead of focusing on the cisternal partitioning of enzymes, we identify the chemical reactions that occur with high propensity in each cisternae. For this we define an effective reaction rate $R_{e f f} (j, k)$ for $P c_{k} \to P c_{k + 1}$ in the jth cisterna as

R_{e f f} (j, k) = \sum_{α = 1}^{N_{E}} R_{α}^{(j)} P^{(j)} (k, α) .

According to our model presented in the section ‘Chemical reaction and transport network in cisternae’, 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 5ai shows the heat map of the effective reaction rates in each cisterna for the optimal that minimizes the normalized KL distance to the 20-GMM target distribution $T^{(20)}$ (see Figure 5aii). The optimized glycan profile displayed in Figure 5aiii is very close to the target. An interesting observation from Figure 5ai is that the same reaction can occur in multiple cisternae.

Figure 5

Download asset Open asset

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-Gaussian mixture model (GMM) target distribution. Here, $N_{E} = 5$ , $N_{C} = 7$ , normalized $D (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.

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 Figure 5bi, we have exchanged the enzymes of the fourth and second cisterna. The glycan profile after enzyme partitioning (see Figure 5biii) is now completely altered (compare Figure 5bii with Figure 5biii). Thus, one can generate different glycan profiles by repartitioning enzymes amongst the same number of cisternae (Jaiman and Thattai, 2018).

Geometry of the fidelity landscape

Here we show that the optimum solution is not unique, rather it is highly degenerate, with several equally good optimum solutions. Thus the multidimensional fidelity landscape in $R$ , $μ$ , $L$ , and $σ$ is typically rugged. We analyse the geometry of this fitness landscape by doing a local Hessian analysis about the optimal solutions.

Degeneracies in the synthesis model

The synthesis model is highly degenerate, in the sense that many combinations of parameters give rise to the same glycan profile. This makes the optimization non-convex as there are many equally good minima. These degeneracies are both discrete and continuous. The continuous degeneracies correspond to regions in reaction rate ( $R$ )-transport rate ( $μ$ ) space moving along which does not change the concentration profile. The discrete degeneracies are disconnected regions in the parameter space which correspond to the same glycan profile. The number of discrete degeneracies increases exponentially with increase in ( $N_{E}, N_{C}$ ). We also find that the fraction of initial conditions converging to a solution close to the global minima increases on increasing ( $N_{E}, N_{C}$ ). Technical details of these issues are discussed in Appendix 8.

Stiff and sloppy directions

We analyse the change in fidelity on small perturbations in $R$ , $μ$ , $L$ , and $σ$ around the optimal solution. This allows us to determine where the cell needs to develop a tighter control mechanism (stiff directions) and where it has more leeway around the optimal values (sloppy directions). We do this by analysing the eigenvalues and eigenvectors of the Hessian around the optimal point (details in Appendix 9). We find that small perturbations around the optimal values in $σ$ change the glycan profile a lot more compared to perturbations in the other parameters and this stiffness in $σ$ generally decreases on increasing $N_{E}, N_{C}$ (Figure 6a–c). Small perturbations in $μ$ and some $L$ directions around the optimum also significantly alter the glycan profile and the stiffness increases on increasing $N_{C}, N_{E}$ , eventually becoming comparable to $σ$ . The glycan profile is robust to perturbations in most $R$ and some $L$ directions (Figure 6b). The total average stiffness of the optimization parameters, defined by the mean of all eigenvalues of the hessian, decreases on increasing $N_{E}, N_{C}$ (Figure 6c).

Figure 6

Download asset Open asset

Stiff and sloppy directions in the optimization parameters.

(a) Eigenvectors of the Hessian matrix $\frac{\partial^{2}}{\partial X_{i} \partial X_{j}} F |_{X_{\min}}$ for $(N_{E}, N_{C}) = (4, 4)$ . The x-axis indexes the $N_{C} + 2 N_{E} N_{C} + 1 = 37$ eigenvectors, the y-axis indexes the $N_{C} + 2 N_{E} N_{C} + 1$ *components* of the eigenvectors, and the greyscale denotes the absolute value of the component in the range $[0, 1]$ . The components are grouped according to ( $μ, R, L, σ$ ), and the eigenvectors are ordered according to the most dominant component in the eigenvector ( $μ$ , orange; $R$ , blue; $L$ , green; $σ$ , purple). There is some mixing of the different components ( $R$ and $μ$ or $σ$ and $μ$ ) but this is usually small. (b) The distribution of eigenvalues $λ_{i}$ of the Hessian matrix $\frac{\partial^{2}}{\partial X_{i} \partial X_{j}} F |_{X_{\min}}$ . Each stripe represents an eigenvalue, and the location of the stripe on the x-axis represents whether the dominant component of the associated eigenvector belongs to $μ$ , $R$ , $L$ , or $σ$ direction. (c) The average stiffness along $μ$ , $R$ , $L$ , or $σ$ directions, defined by the log of the average of eigenvalues corresponding to the eigenvectors in the respective group, as a function of $N_{C}$ for fixed $N_{E} = 4$ . (d) Total average stiffness $⟨ λ ⟩ = \log (\frac{\sum λ_{i}}{N_{C} + 2 N_{E} N_{C} + 1})$ as a function of $N_{E}, N_{C}$ .

Implications for robustness to parametric noise

Since the synthesized glycan distribution displayed by the cell marks its identity, it must be robust to noise intrinsic to the synthesis machinery. The degeneracy of solutions and sloppy directions in the fidelity landscape makes the glycan distribution robust to intrinsic noise in the synthesis and cell-to-cell variations in the kinetic parameters. We find that the number of degeneracies increases on increasing ( $N_{E}, N_{C}$ ), and the average stiffness of the optimized parameters decreases on increasing ( $N_{E}, N_{C}$ ), making the synthesis more robust to parameter fluctuations. Further, while the parameter space is high dimensional, the dimension of controllable parameters (measured by the stiff directions) is low dimensional. We find this dimensional reduction a compelling idea which we will take up later.

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 the diversity of glycans (Varki, 2011; Dennis et al., 2009). 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 Brunner’s gland cells secrete mucous that contains heavily O-glycosylated mucin as its major component (Van Halbeek et al., 1983). 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 (McFarlane et al., 2014; Koch et al., 2015; O’Neill et al., 2004; Hayashi and Kaida, 2011; Kumar et al., 2011; Gow and Hube, 2012; Atmodjo et al., 2013; Free, 2013; Pauly et al., 2013; Burton and Fincher, 2014). These organisms possess enlarged Golgi complexes with multiple cisternae per stack (Becker and Melkonian, 1996; Mironov et al., 2017; Donohoe et al., 2007; Mogelsvang et al., 2003; Ladinsky et al., 2002).

We define diversity as the total number of glycan species produced above a specified threshold abundance $c_{t h}$ . 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_{t h} = 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 (see Appendix 10—figure 1).

We use the sigmoid function ${(1 + e^{- x / τ})}^{- 1}$ as a differentiable approximation to the Heaviside function $Θ (x)$ and define the following optimization to maximize diversity for a given set of parameter values, $N_{E}, N_{C}, σ$ :

Diversity (σ, N_{C}, N_{E}) := \begin{aligned} max_{μ, R, L} & \sum_{i = 1}^{N_{s}} (1 + e^{- N_{s} (c_{i} - c_{t h})})^{- 1} \\ s.t. & R_{min} \leq R_{α}^{(j)} \leq R_{max}, \\ μ_{min} \leq μ^{(j)} \leq μ_{max}, \end{aligned}

where, as before, $(μ_{\max}, μ_{\min}) = (1, 0.01)$ /min, and $(R_{\max}, R_{\min}) = (20, 0.018)$ /min, and $c_{t h} = 1 / N_{s}$ is the threshold. See Appendix 2 for details on the parameter estimation.

The results displayed in Figure 7a 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 correlational at best since there are many ways to achieve high glycan diversity – notably by increasing the number of enzymes.

Figure 7

Download asset Open asset

Strategies for achieving high glycan diversity.

Diversity versus $N_{C}$ and transport rate $μ$ at various values of specificity $σ$ for fixed $N_{E} = 3$ . (a) Diversity vs. $N_{C}$ at optimal transport rate $μ$ . Diversity initially increases with $N_{C}$ , but eventually levels off. The levelling off starts at a higher $N_{C}$ 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_{\min}^{- 1}$ ) at various value of $σ$ , for fixed $N_{C} = 4$ and $N_{E} = 10$ .

On the other hand, one of the goals of glycoengineering is to produce a particular glycan profile with low heterogeneity (Fisher et al., 2019; Jaiman and Thattai, 2018). 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 the 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, which is consistent with the results in Fisher et al., 2019. To get a sharper peak, as advocated for instance by Jaiman and Thattai, 2018, one might need to increase the number of high-specificity enzymes $N_{E}$ further.

Discussion

The precision of the stereochemistry and enzymatic kinetics of these N-glycosylation reactions (Varki, 2009) has inspired a number of mathematical models (Umaña and Bailey, 1997; Krambeck et al., 2009; Krambeck and Betenbaugh, 2005) 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 (Umaña and Bailey, 1997; Krambeck and Betenbaugh, 2005; Krambeck et al., 2009) 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 (Fisher et al., 2019; Fisher and Ungar, 2016) 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 can help predict glycosylation patterns and direct glycoengineering projects to optimize glycoform distributions.

Our focus is different. We are interested in the role of glycans as a marker or molecular code of cell identity (Gabius, 2018; Varki, 2017; Pothukuchi et al., 2019), and in particular, understanding enzymatic and transport processes located in the secretory apparatus of the cell that ensure that this code is generated with high fidelity. To do this, we have had to develop a new formal apparatus that allows us to address these questions and discuss trade-offs between competing drives. Since our analysis draws on many diverse fields, we provide a short summary of the assumptions, methods, and results of the article before discussing the implications of our work.

The glycan profile on the cell surface is a marker of cell-type identity (Varki, 2009; Gabius, 2018; Varki, 2017; Pothukuchi et al., 2019). We define the complexity of a glycan profile to be the minimum number of GMM components required to approximate the profile to within the noise floor. We show that with this definition of complexity more complex organisms correlate with higher complexity glycan profiles. We use this to analyse the complexity of the glycan profiles of planaria, hydra, and mammalian cells (Drickamer and Taylor, 1998).
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 (Varki, 2017; Varki, 1998; Pothukuchi et al., 2019). We have proposed a general model for chemical synthesis and transport that, in principle, allows us to compute the synthesized glycan distribution at the cell surface as a function of the enzymes $N_{E}$ , reaction rates $R$ , enzyme configurations $L$ , specificity of enzymes $σ$ , number of cisternae $N_{C}$ , and transport rates $μ$ . However the large dimension of the search space makes this optimization intractable. We thus use a simplified synthesis model with fewer parameters; while our quantitative results are based on this simplified model, we believe that at a qualitative level our results have more general validity.
We define the fidelity of a synthesis mechanism as the minimum normalized KL divergence (Cover and Thomas, 2012; MacKay, 2003) between synthesized glycan distribution on the cell surface and a ‘target’ profile.
The results of the optimization over rates $R$ and enzyme configurations $L$ for a given value of $(N_{C}, N_{E}, σ)$ and a target distribution $c^{*}$ of given complexity are given in Figure 3 and Figure 4. Here, we highlight some qualitative consequences of the model:
1. Keeping the number of enzymes fixed, a more elaborate transport mechanism (via control of $N_{C}$ and $μ$ ) is essential for synthesizing high-complexity target distributions to within a high fidelity, or equivalently, low error (Figure 4a and b). Fewer cisternae cannot be compensated for by optimizing the enzymatic synthesis via control of parameters $R$ , $L$ , and $σ$ . An empirical verification of this would involve a coordinated analysis of the glycan profiles, ultrastructure of Golgi, and the number of glycosylation enzymes across many species.
2. Thus, our study suggests that the requirement that a glycan code of a given complexity be synthesized with sufficiently high fidelity imposes functional control on the Golgi cisternal number. It also provides an argument for the evolutionary requirement of multiple compartments by demonstrating that the fidelity and robustness of the glycan code arising from a chemical synthesis that involves multiple cisternae are higher than the one that involves a single cisterna (keeping everything else fixed) (see Figure 4a and b and Figure 6) 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 Jaiman and Thattai, 2018.
  
  Combining (a) and (b), our study quantitatively shows that constructing a high-fidelity representation of a complex target distribution, such as those observed in real cells, requires a complex Golgi machinery with multiple cisternae, precise enzyme partitioning, and control on enzyme specificity. This definition of fidelity of the glycan code allows us to provide a quantitative argument for the evolutionary requirement of multiple compartments. While it is possible to produce complex glycan distributions in one compartment using a large number of enzymes, such a design would inevitably require a more elaborate genetic cost.
3. Organisms such as plants and algae have a diverse repertoire of glycans that are utilized in a variety of functions (McFarlane et al., 2014; Koch et al., 2015; O’Neill et al., 2004; Hayashi and Kaida, 2011; Kumar et al., 2011; Gow and Hube, 2012; Atmodjo et al., 2013; Free, 2013; Pauly et al., 2013; Burton and Fincher, 2014). Our study shows that it is optimal to use low-specificity enzymes to synthesize target distributions with high diversity (Figure 7). 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}$ .
4. 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 Figure 4d, this optimal enzyme specificity can be very high for highly complex target distributions. Such high specificity can lower fitness when the environment, and hence the target glycan distribution, fluctuates rapidly, and the synthesis parameters cannot change rapidly enough to track the environment (Nam et al., 2012; Peracchi, 2018). This compromise, between robustness to a changing environment and high fidelity in synthesizing high-complexity glycan profiles, is achievable by sloppy enzymes coupled with error-correcting mechanisms (Nam et al., 2012; Peracchi, 2018). However, sloppy enzymes create ‘wrong’ glycans, and therefore, ex-post error-correcting mechanisms must be in place to correct synthesis errors to ensure high fidelity of the glycan code. A task for the future is to understand the role of intracellular transport in providing non-equilibrium proofreading mechanisms to reduce such coding errors, and its optimal adaptive strategies and plasticity in a time-varying environment.
  
  Combining (c) and (d), we find that 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. Sloppy or promiscuous enzymes bring in the potential for evolvability (Kirschner and Gerhart, 2008), and sloppiness allows the system to be stable to random mutations in proteins or variations in the target distribution.
5. The model solution is degenerate, in the sense that there are many equally good global minimas. These degeneracies are both continuous and discrete. The continuous degeneracies correspond to regions in the reaction rate – transport rate space, moving along which will not change the concentration profile, thus ensuring robustness to internal noise. This suggests that the distribution is robust to slight cell-to-cell variations in these kinetic parameters.
6. Our model implies that close to a local minima the inter-cisternal transport rate $μ$ and the specificity of the enzymes $σ$ are stiff directions, that is, the cell should exercise tighter control on $μ$ and $σ$ as compared to the other parameters. The reaction rates close to the local minima are sloppy directions, and moving along these directions does not change the glycan profile much.
Taken together, our quantitative analysis of the trade-offs has deep implications for non-equilibrium self-assembly of the Golgi cisternae, and suggests that the non-equilibrium control of cisternal number must involve a coupling of non-equilibrium self-assembly of cisternae with enzymatic chemical reaction kinetics (Glick and Malhotra, 1998).

Admittedly the chemical network that we have considered here is much simpler than the chemical network associated with the 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 (Alberts, 2002; Varki, 2009; Pothukuchi et al., 2019). Our task has been to get at a qualitative understanding using quantitative methods and thereby to arrive at general principles. We believe our analysis is generalizable 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 high-fidelity synthesis of a complex glycan code. 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 (Sengupta and Linstedt, 2011). An important implication is that a description of the non-equilibrium self-assembly of a fixed number of Golgi cisternae must combine the dynamics of chemical processing and membrane dynamics involving fission, fusion, and transport (Sengupta and Linstedt, 2011; Sachdeva et al., 2016; Sens and Rao, 2013). We believe that this is a promising direction for future research.

Appendix 1

Constructing target distributions for glycans of a given cell type

The distribution of the glycans on the cell surface is obtained via mass spectrometry. The x-axis of mass spectroscopy (MS) graphs is mass/charge of the ionized sample molecules, and the y-axis is relative intensity corresponding to each mass/charge value, taking the highest intensity as 100%. This relative intensity roughly correlates with the relative abundances of the molecules in the sample.

The raw MS data is noisy and cannot be directly used as the target distribution in our optimization problem. There are three major sources of noise in the MS data (Du et al., 2008): chemical noise in the sample, the Poisson noise associated with detecting discrete events, and the Nyquist–Johnson noise associated with any charge system. We propose a simple model that accounts for the chemical noise and the Poisson sampling noise. Using this noise model and the available MS data, we generate parametric bootstrap samples of glycan measurements and fit a GMM on this sample to approximate the glycan distribution. This GMM probability distribution is used as the target distribution in our numerical experiments.

The MS data obtained from Cummings and Crocker, 2020; Subramanian et al., 2018; Sahadevan et al., 2014 had mass ranging between 500 and 5000 Da with intensity reported at every 0.0153 Da. We first bin this MS data into 180 bins and take the maximum value within each bin as the value of intensity for that bin. Appendix 1—figure 1 plots the raw MS data and the binned distribution. Next, we describe the parametric bootstrap model that we used to generate the glycan data. Let ${\bar{I}}_{k}$ represent the relative intensity of the kth bin in the binned MS graph. We generate a sample population of glycans using the MS data in the following way:

Poisson sampling noise: The MS data does not have absolute count information. We assume an arbitrary maximum count $I_{\max}$ and define the intensity $I_{k} = I_{\max} {\bar{I}}_{k}$ . The plots in Appendix 1—figure 2a show that the results are not sensitive to the specific value of $I_{\max}$ .
Chemical noise: The sample used for MS analysis also contains small amounts of molecules that are not glycans. These appear as very small peaks in the MS data. We assume that the probability p_k that the peak at index $k$ corresponds to a glycan is given by
$p_{k} = 1 - e^{- \frac{I_{k}}{I_{m a x}}} = 1 - e^{- {\bar{I}}_{k}}$
which adequately suppresses this chemical noise.
Bootstrapped glycan data: The count n_k at the glycan index $k$ is distributed according to the following distribution:
$n_{k} = {\begin{matrix} 0 & (1 - p_{k}) & n = 0 \\ n & p_{k} e^{- I_{k}} \frac{{(I_{k})}^{n}}{n!} & n \geq 1 . \end{matrix}$
We assume that the MS data was generated from N different cells. Thus, the total count at glycan index k is given by the sum of N i.i.d. samples distributed according to the distribution above. In Appendix 1—figure 2b, we show that results are insensitive to N. We normalize the count distribution by the total number of counts across all the bins to obtain the bootstrapped probability mass function PT.

Appendix 1—figure 1

Download asset Open asset

The binned mass spectroscopy (MS) data (blue) approximates the raw MS data (red) very well.

We use this binned data for Gaussian mixture model (GMM) approximation of the MS data.

Appendix 1—figure 2

Download asset Open asset

Log likelihood vs. number of components ( $N$ ) in the Gaussian mixture model (GMM).

We see that the log likelihood saturates at around, $N = 20$ , thus 20-GMM is a very good representation of the mass spectroscopy (MS) data from *human* T-cells. The different symbols are for (a) different values of the maximum intensity $I_{m a x} = 50, 100, 200$ and (b) different values of the number of i.i.d. samples, $N_{i} = 500, 1000, 2000$ showing the insensitivity of the log likelihood to the value of $I_{m a x}$ and. $N_{i}$

The bootstrapped distribution $p_{T}$ is noisy, and hence cannot be used directly as the target distribution. We use a GMM-based approach to denoise the raw data. The advantage of using a GMM-based approach is that it creates an easily interpretable hierarchy of increasingly more detailed distributions to approximate the mass spectrometry profile. We define the complexity of a mass spectrometry profile as the minimum number of components (individual Gaussians) in the GMM model required to approximate it. The details of the GMM calculations are as follows. We fix the number of components $m$ . We want to approximate the bootstrapped probability $p_{T}$ by the $m$ -component mixture of Gaussian distributions $p_{G M M} (θ) = \sum_{i = 1}^{m} w_{i} N_{η_{i}, Δ_{i}}$ , where $N_{η_{i}, Δ_{i}}$ denotes the Gaussian distribution with mean $η_{i}$ and variance $Δ_{i}$ , $w_{i} \geq 0$ and $\sum_{i = 1}^{m} w_{i} = 1$ , the parameter vector $θ = (w, η, Δ)$ . We compute the optimal $m$ -component GMM approximation by minimizing the KL divergence $D (p_{T} | | p_{G M M} (θ))$ as a function of parameter vector $θ$ . Since

\begin{array}{ll} D (p_{T} | | p_{G M M} (θ)) & := \sum_{k = 1}^{N_{s}} p_{T} (k) \log (\frac{p_{T} (k)}{\sum_{i}^{m} w_{i} N_{η_{i}, Δ_{i}} (k)}) \\ = \sum_{k = 1}^{N_{s}} p_{T} (k) \log p_{T} (k) - \sum_{k = 1}^{N_{s}} p_{T} (k) \log (\sum_{i}^{m} w_{i} N_{η_{i}, Δ_{i}} (k)), \end{array}

the optimization problem $\min_{θ} D (p_{T} | | p_{G M M} (θ))$ is equivalent to

max_{θ} g (θ) := \sum_{k = 1}^{N_{S}} p_{T} (k) \log (\sum_{i = 1}^{m} w_{i} N_{η_{i}, Δ_{i}} (k))

This is a non-convex optimization. We use an expectation-maximization (EM)-based iterative heuristic to compute a local maximum. Let $θ^{(t)}$ denote the current value of the parameters. For each component $i = 1, \dots, m$ , and index $k = 1, \dots, N_{s}$ , define

z_{i}^{(t)} (k) = \frac{w_{i}^{(t)} N_{η_{i}^{(t)}, Δ_{i}^{(t)}} (k)}{\sum_{j = 1}^{m} w_{j}^{(t)} N_{η_{j}^{(t)}, Δ_{j}^{(t)}} (k)} .

Then $z_{i}^{(t)} (k) \geq 0$ , and $\sum_{i = 1}^{m} z_{i}^{(t)} (k) = 1$ . We interpret $z_{i}^{(t)} (k)$ as the probability that the count in bin $k$ came from component $i$ . Define

Q (θ, θ^{(t)}) = \sum_{k = 1}^{N_{S}} \sum_{i = 1}^{m} p_{T} (k) z_{i}^{(t)} (k) \log (\frac{w_{i} N_{η_{i}, Δ_{i}} (k)}{z_{i}^{(t)} (k)})

Then, we have that

Q (\hat{θ}, \hat{θ}) = \sum_{k = 1}^{N_{S}} \sum_{i = 1}^{m} p_{T} (k) {\hat{z}}_{i} (k) \log (\sum_{i = 1}^{m} w_{i} N_{η_{i}, Δ_{i}} (k)) = \sum_{k = 1}^{N_{S}} p_{T} (k) \log (\sum_{i = 1}^{m} w_{i} N_{η_{i}, Δ_{i}} (k)) = g (\hat{θ}),

and

g (θ) = \sum_{k = 1}^{N_{S}} p_{T} (k) \log (\sum_{i = 1}^{m} \frac{w_{i} N_{η_{i}, Δ_{i}} (k)}{z_{i}^{(t)} (k)} z_{i}^{(t)} (k)) \geq \sum_{k = 1}^{N_{S}} \sum_{i = 1}^{m} p_{T} (k) z_{i}^{(t)} (k) \log (\frac{w_{i} N_{η_{i}, Δ_{i}} (k)}{z_{i}^{(t)} (k)}) = Q (θ, θ^{(t)}) .

Define

θ^{(t + 1)} = arg max_{θ} Q (θ, \hat{θ})

Then, we have that

g (θ^{(t + 1)}) \geq Q (θ^{(t + 1)}, θ^{(t)}) \geq Q (θ^{(t)}, θ^{(t)}) = g (θ^{(t)}) .

Therefore, the iterative algorithm in (Equation 11) generates a sequence ${θ^{(t)} : t \geq 1}$ with non-decreasing values of $g$ , and the sequence converges to a local maximum. Next, we show that the optimization in (8) can be computed efficiently.

$w$ -update
$w^{(t + 1)} = arg max_{w} \sum_{k = 1}^{N_{S}} \sum_{i = 1}^{m} p_{T} (k) z_{i}^{(t)} (k) \log (w_{i}) ⟹ w_{i}^{(t + 1)} = \frac{\sum_{k = 1}^{N_{S}} z_{i}^{(t)} (k) p_{T} (k)}{\sum_{i = 1}^{m} \sum_{k = 1}^{N_{S}} z_{i}^{(t)} (k) p_{T} (k)}$
$η$ -update
$η_{i}^{(t + 1)} = arg min_{η} \sum_{k = 1}^{N_{S}} p_{T} (k) {\hat{z}}_{i} (k) | k - η_{i} |^{2} ⟹ η_{i}^{(t + 1)} = \frac{\sum_{k = 1}^{N_{s}} z_{i}^{(t)} (k) k}{\sum_{k = 1}^{N_{s}} z_{i}^{(t)} (k)} .$
$Δ$ -update
$\begin{array}{ll} Δ_{i}^{(t + 1)} & = {argmax}_{Δ \geq Δ_{c u t}} {- \frac{\sum_{k = 1}^{N_{s}} p_{T} (k) z_{i}^{(t)} (k) | k - η_{i}^{t + 1} |^{2}}{2 Δ} - \log (Δ)} \\ = max (\sqrt{\frac{\sum_{k = 1}^{N_{S}} p_{T} (k) z_{i}^{(t)} (k) | k - η_{i}^{(t + 1)} |^{2}}{\sum_{k = 1}^{N_{s}} p_{T} (k) z_{i}^{(t)} (k)}}, Δ_{c u t}), \end{array}$

where $Δ_{c u t}$ is the minimum allowed width of the Gaussians, in our case $Δ_{c u t} = 1$ since glycan index, $k \in {1, 2, \dots N_{S}}$ , takes integer values with spacing 1.

Since this is a heuristic algorithm for a non-convex optimization, we performed several initializations of the algorithm to identify the best local maximum.

The number of components $m$ in a GMM is a free parameter. The KL divergence between the true and GMM approximated ( $D (p_{T} | | p_{G M M})$ ), shown in Figure 1, saturates at some value of the number of components, and adding components beyond this only increases model complexity without increasing the quality of approximation. We define the complexity of a mass spectrometry data by the number of components at which saturation is reached. We compare the complexity of glycan profiles of hydra, planaria, and humans. The numbers of cell types in hydra, planaria, and humans are around 41 (Siebert et al., 2019), 44 (Fincher et al., 2018), and 103 (Han et al., 2020), respectively, based on transcriptome analysis (these are lower bounds based on the main cell types, and especially for planaria and hydra, are subject to constant revision). Our analysis of the MSMS data of these organisms suggests that organisms with fewer cell types have less complex glycan distribution.

Appendix 2

Parameter estimation

The typical transport time of glycoproteins across the Golgi complex is estimated to be in the range 15–20 mins (Umaña and Bailey, 1997), which corresponds to the transport rate $μ = 0.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⁶ cell 24 hr (Umaña and Bailey, 1997; Krambeck et al., 2009). For our calculation, we set $q = 387.30$ pmol/10⁶ cells 24 hr = 0.27 pmol/10⁶ cells min, where 387.30 is the geometric mean of 100 and 1500. We set the range for the enzymatic rate $R$ to be

R_{min} = min_{α} {\frac{V_{max}^{(α)} / ν}{K_{M}^{(α)} + \frac{1}{ν} \frac{q}{μ}}} \leq R \leq R_{max} = max_{α} {\frac{V_{max}^{(α)} / ν}{K_{M}^{(α)}}} .

where $K_{M}^{(α)}$ and $V_{\max}^{(α)}$ denote the Michaelis constants and $V_{\max}$ of the αth enzyme. The conversion from 1 pmol/10⁶ cells to concentration can be obtained by taking cisternal volume ( $ν$ ) to be 2.5 µm³ (Umaña and Bailey, 1997; Krambeck et al., 2009). This gives

1 pmoles/ 10^{6} cells = \frac{10^{- 12} moles}{10^{6} \times 2.5 \times 10^{- 18} \times 10^{3} litre} = 400 μ M .

In Appendix 2—table 1, we report the parameters for the eight enzymes taken from Table 3 in Umaña and Bailey, 1997. From these parameters, it follows that

\begin{array}{ll} R_{min} & = min_{α} {\frac{V_{max}^{(α)} / ν}{K_{M}^{(α)} + \frac{1}{ν} \frac{q}{μ}}} \\ = \frac{V_{max}^{(7)} / ν}{K_{M}^{(7)} + \frac{1}{ν} \frac{q}{μ}} = \frac{.16 \times 400 μ M / min}{3400 μ M + 149.4 μ M} = 0.018 {min}^{- 1} \end{array}

\begin{array}{ll} R_{max} & = max_{α} {\frac{V_{max}^{(α)} / ν}{K_{M}^{(α)}}} \\ = \frac{V_{max}^{(1)} / ν}{K_{M}^{(1)}} = \frac{5 \times 400 μ M / min}{100 μ M} = 20 {min}^{- 1} \end{array}

Appendix 2—table 1

Enzyme parameters taken from Table 3 in Umaña and Bailey, 1997 that we use to calculate the bounds on the reaction rate $R$ .

Here, $K_{M}^{(α)}$ and $V_{\max}^{(α)}$ denote the Michaelis constant and $V_{\max}$ of the αth enzyme.

$α$	$K_{M}^{(α)}$ (µmol)	$V_{\max}^{(α)}$ (pmol/10⁶ cell-min)
$α$	$K_{M}^{(α)}$ (µmol)	$V_{\max}^{(α)}$ (pmol/10⁶ cell-min)	1	100	5
2	260	7.5
3	200	5
4	100	5
5	190	2.33
6	130	.16
7	3400	.16
8	4000	9.66

Appendix 3

Kinetics of sequential chemical reactions and transport

On including the rates of injection ( $q$ ) and inter-cisternal transport $μ^{(j)}$ from the cisterna $j$ to $j + 1$ into the chemical reaction kinetics, the change in substrate concentrations $c_{k}^{(j)}$ with time is given by

\begin{array}{ll} \frac{d c_{1}^{(1)}}{d t} = & q - \sum_{α = 1}^{N_{E}} \frac{V (1, 1, α) P^{(1)} (1, α) c_{1}^{(1)}}{M (1, 1, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(1)} (k^{'}, α) c_{k^{'}}^{(1)}}{M (1, k^{'}, α)})} - μ^{(1)} c_{1}^{(1)} \\ \frac{d c_{k}^{(1)}}{d t} = & \sum_{α = 1}^{N_{E}} \frac{V (1, k - 1, α) P^{(1)} (k - 1, α) c_{k - 1}^{(1)}}{M (1, k - 1, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(1)} (k^{'}, α) c_{k^{'}}^{(1)}}{M (1, k^{'}, α)})} \\ - \sum_{α = 1}^{N_{E}} \frac{V (1, k, α) P^{(1)} (k, α) c_{k}^{(1)}}{M (1, k, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(1)} (k^{'}, α) c_{k^{'}}^{(1)}}{M (1, k^{'}, α)})} - μ^{(1)} c_{k}^{(1)} \\ \frac{d c_{N_{s}}^{(1)}}{d t} = & \sum_{α = 1}^{N_{E}} \frac{V (1, N_{s} - 1, α) P^{(1)} (N_{s} - 1, α) c_{N_{s} - 1}^{(1)}}{M (1, N_{s} - 1, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(1)} (k^{'}, α) c_{k^{'}}^{(1)}}{M (1, k^{'}, α)})} - μ^{(1)} c_{N_{s}}^{(1)} \end{array}

for cisterna-1, and

\begin{array}{ll} \frac{d c_{1}^{(j)}}{d t} = & μ^{(j - 1)} c_{1}^{(j - 1)} - \sum_{α = 1}^{N_{E}} \frac{V (j, 1, α) P^{(j)} (1, α) c_{1}^{(j)}}{M (j, 1, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(j)} (k^{'}, α) c_{k^{'}}^{(j)}}{M (j, k^{'}, α)})} - μ^{(j)} c_{1}^{(j)} \\ \frac{d c_{k}^{(j)}}{d t} = & μ^{(j - 1)} c_{k}^{(j - 1)} + \sum_{α = 1}^{N_{E}} \frac{V (j, k - 1, α) P^{(j)} (k - 1, α) c_{k - 1}^{(j)}}{M (j, k - 1, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(j)} (k^{'}, α) c_{k^{'}}^{(j)}}{M (j, k^{'}, α)})} \\ - \sum_{α = 1}^{N_{E}} \frac{V (j, k, α) P^{(j)} (k, α) c_{k}^{(j)}}{M (j, k, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(j)} (k^{'}, α) c_{k^{'}}^{(j)}}{M (j, k^{'}, α)})} - μ^{(j)} c_{k}^{(j)} \\ \frac{d c_{N_{s}}^{(j)}}{d t} = & μ^{(j - 1)} c_{N_{s}}^{(j - 1)} + \sum_{α = 1}^{N_{E}} \frac{V (j, N_{s} - 1, α) P^{(j)} (N_{s} - 1, α) c_{N_{s} - 1}^{(j)}}{M (j, N_{s} - 1, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(j)} (k^{'}, α) c_{k^{'}}^{(j)}}{M (j, k^{'}, α)})} - μ^{(j)} c_{N_{s}}^{(j)} \end{array}

for cisternae $j = 2, 3, \dots, N_{C}$ . These set of dynamical equations (13)-(14), with initial conditions, can be solved to obtain the concentration $c_{k}^{(j)} (t)$ for $t \geq 0$ . Equations (13)-(14) automatically obey the conservation law for the protein concentration ( $p$ ), that is, the total protein concentration $p^{(j)} = \sum_{k^{'} = 1}^{N_{s}} c_{k^{'}}^{(j)}$ in the jth cisterna automatically satisfies

\begin{array}{ll} \frac{d p^{(1)}}{d t} = q - μ^{(1)} p^{(1)} \\ \frac{d p^{(j)}}{d t} = μ^{(j - 1)} p^{(j - 1)} - μ^{(j)} p^{(j)} \end{array}

for $j = 2, 3, \dots N_{C}$ .

At steady state, the left-hand side of equations (13)-(14) is set to zero, which, after rescaling the kinetic parameters in terms of the injection rate $q$ , that is, $V (j, k, α) = V (j, k, α) / q$ and $μ^{(j)} = μ^{(j)} / q$ , gives the following recursion relations for the steady-state concentrations of the glycans in each cisterna. In the first cisterna,

\begin{array}{ll} c_{1}^{(1)} = \frac{1}{μ^{(1)} + \sum_{α = 1}^{N_{E}} \frac{V (1, 1, α) P^{(1)} (1, α) c_{1}^{(1)}}{M (1, 1, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(1)} (k^{'}, α) c_{k^{'}}^{(1)}}{M (1, k^{'}, α)})}} \\ c_{k}^{(1)} = \frac{\sum_{α = 1}^{N_{E}} \frac{V (1, k - 1, α) P^{(1)} (k - 1, α) c_{k - 1}^{(1)}}{M (1, k - 1, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(1)} (k^{'}, α) c_{k^{'}}^{(1)}}{M (1, k^{'}, α)})}}{μ^{(1)} + \sum_{α = 1}^{N_{E}} \frac{V (1, k, α) P^{(1)} (k, α) c_{k}^{(1)}}{M (1, k, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(1)} (k^{'}, α) c_{k^{'}}^{(1)}}{M (1, k^{'}, α)})}} \\ c_{N_{s}}^{(1)} = \frac{\sum_{α = 1}^{N_{E}} \frac{V (1, N_{s} - 1, α) P^{(1)} (N_{s} - 1, α) c_{N_{s} - 1}^{(1)}}{M (1, N_{s} - 1, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(1)} (k^{'}, α) c_{k^{'}}^{(1)}}{M (1, k^{'}, α)})}}{μ^{(1)}} \end{array}

and in cisternae $j \geq 2$ ,

\begin{array}{ll} c_{1}^{(j)} = \frac{μ^{(j - 1)} c_{1}^{(j - 1)}}{μ^{(j)} + \sum_{α = 1}^{N_{E}} \frac{V (j, 1, α) P^{(j)} (1, α) c_{1}^{(j)}}{M (j, 1, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(j)} (k^{'}, α) c_{k^{'}}^{(j)}}{M (j, k^{'}, α)})}} \\ c_{k}^{(j)} = \frac{μ^{(j - 1)} c_{k}^{(j - 1)} + \sum_{α = 1}^{N_{E}} \frac{V (j, k - 1, α) P^{(j)} (k - 1, α) c_{k - 1}^{(j)}}{M (j, k - 1, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(j)} (k^{'}, α) c_{k^{'}}^{(j)}}{M (j, k^{'}, α)})}}{μ^{(j)} + \sum_{α = 1}^{N_{E}} \frac{V (j, k, α) P^{(j)} (k, α) c_{k}^{(j)}}{M (j, k, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(j)} (k^{'}, α) c_{k^{'}}^{(j)}}{M (j, k^{'}, α)})}} \\ c_{N_{s}}^{(j)} = \frac{μ^{(j - 1)} c_{N_{s}}^{(j - 1)} + \sum_{α = 1}^{N_{E}} \frac{V (j, N_{s} - 1, α) P^{(j)} (N_{s} - 1, α) c_{N_{s} - 1}^{(j)}}{M (j, N_{s} - 1, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(j)} (k^{'}, α) c_{k^{'}}^{(j)}}{M (j, k^{'}, α)})}}{μ^{(j)}} \end{array}

Equations (15)-(16) automatically imply that the total steady-state glycan concentration in each cisterna $j = 1, \dots, N_{c}$ is given by

\sum_{k = 1}^{N_{s}} c_{k}^{(j)} = \frac{1}{μ^{(j)}} .

Appendix 4

Reformulation of Optimization A

Define a new set of parameters,

R (j, k, α) = \sum_{α = 1}^{N_{E}} \frac{V (j, k, α)}{M (j, k, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(j)} (k^{'}, α) c_{k^{'}}^{(j)}}{M (j, k^{'}, α)})}

where $c$ denotes the steady-state glycan concentration corresponding to a specific $(μ, M, V, L)$ . Define $v$ by the following set of linear equations:

\begin{array}{ll} v_{1}^{(1)} = & \frac{1}{μ^{(1)} + \sum_{α = 1}^{N_{E}} R (1, 1, α) P^{(1)} (1, α)} \\ v_{k}^{(1)} = & \frac{v_{k - 1}^{(1)} \sum_{α = 1}^{N_{E}} R (1, k - 1, α) P^{(1)} (k - 1, α)}{μ^{(1)} + \sum_{α = 1}^{N_{E}} R (1, k, α) P^{(1)} (k, α)} \\ v_{N_{s}}^{(1)} = & \frac{v_{N_{s} - 1}^{(1)} \sum_{α = 1}^{N_{E}} R (1, N_{s} - 1, α) P^{(1)} (N_{s} - 1, α)}{μ^{(1)}} \end{array}

for $j = 1$ , and

\begin{array}{ll} v_{1}^{(j)} = & \frac{v_{1}^{(j - 1)} μ^{(j - 1)}}{μ^{(j)} + \sum_{α = 1}^{N_{E}} R (j, 1, α) P^{(j)} (1, α)} \\ v_{k}^{(j)} = & \frac{v_{k}^{(j - 1)} μ^{(j - 1)}}{μ^{(j)} + \sum_{α = 1}^{N_{E}} R (j, k, α) P^{(j)} (k, α)} \\ + \frac{v_{k - 1}^{(j)} \sum_{α = 1}^{N_{E}} R (j, k - 1, α) P^{(j)} (k - 1, α)}{μ^{(j)} + \sum_{α = 1}^{N_{E}} R (j, k, α) P^{(j)} (k, α)} \\ v_{N_{s}}^{(j)} = & \frac{v_{N_{s}}^{(j - 1)} \sum_{α = 1}^{N_{E}} R (j, N_{s} - 1, α) P^{(j)} (N_{s} - 1, α)}{μ^{(j)}} + \frac{v_{N_{s}}^{(j - 1)} μ^{(j - 1)}}{μ^{(j)}} \end{array}

for $j = 2, \dots, N_{C}$ . Then, by the definition of $R$ in (17), it trivially follows that the steady-state concentration $c$ corresponding to $(μ, M, V, L)$ is a solution for (18)–(19).

Next, we show that for $v$ obtained from (18)–(19) for any parameter $(μ, R, L)$ , there exists parameter $(μ, M, V, L)$ such that (15)–(16) are automatically satisfied when we set $c = v$ , that is, $v$ is the steady-state concentration for $(μ, M, V, L)$ , and vice versa. Let

A = {[c_{k}^{(j)}]_{j, k} : \begin{array}{l} μ^{(j)} \geq 0, M (j, k, α) \geq 0, V (j, k, α) \geq 0, l_{α}^{(j)} \geq 0, \\ [c_{k}^{(j)}]_{j k} given by (15) and (16), \end{array}}

and let

B = {[v_{k}^{(j)}]_{j, k} : \begin{array}{l} μ^{(j)} \geq 0, R (j, k, α) \geq 0, l_{α}^{(j)} \geq 0 \\ [v_{k}^{(j)}]_{j, k} given by (18) and (19) \end{array}} .

Then, our task is to show that $A = B$ . Suppose $[c_{k}^{(j)}]_{j, k} \in A$ . Let $[M (j, k, α)]$ , $[V (j, k, α)]$ , and $[l_{α}^{(j)}]$ be the corresponding parameters. Define

R (j, k, α) = \sum_{α = 1}^{N_{E}} \frac{V (j, k, α)}{M (j, k, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(j)} (k^{'}, α) c_{k^{'}}^{(j)}}{M (j, k^{'}, α)})} \geq 0

Then ${[c_{k}^{(j)}]}_{j, k} \in ℬ$ .

Next, suppose ${[v_{k}^{(j)}]}_{j, k} \in ℬ$ . Let $[R (j, k, α)]$ , $[l_{α}^{(j)}]$ denote the corresponding parameters. Since $\sum_{k = 1}^{N_{s}} v_{k}^{(j)} = 1 / μ^{(j)} < \infty$ , it follows that $\sum_{k = 1}^{N_{s}} P^{(j)} (k, α) v_{k}^{(j)} < 1 / μ^{(j)} < \infty$ . Thus, there exist parameters $[M (j, k, α)]$ , $[V (j, k, α)]$ , and $[l_{α}^{(j)}]$ such that

R (j, k, α) = \sum_{α = 1}^{N_{E}} \frac{V (j, k, α)}{M (j, k, α) (1 + \sum_{k^{'} = 1}^{N_{s}} \frac{P^{(j)} (k^{'}, α) v_{k^{'}}^{(j)}}{M (j, k^{'}, α)})}

Therefore, ${[v_{k}^{(j)}]}_{j, k}$ satisfy (15) and (16), that is, $[v_{k}^{(j)}]_{j, k} \in A$ .

Thus, the set of all concentration profiles defined by (18)–(19) as a function of all possible values of the parameters $(μ, R, L)$ is identical to the set defined by (15)–(16) as a function of $(μ, M, V, L)$ . This is a crucial insight since it allows us to search the entire parameter space using (18)–(19), where the concentration is known explicitly in terms of $(μ, R, L)$ .

To pose the new optimization problem, it is convenient to define ${\bar{v}}_{i}^{(j)} = μ^{(j)} v_{i}^{(j)}$ , for $j = 1, \dots, N_{C}$ , and $i = 1, \dots, N_{s}$ . Furthermore, the glycan distribution of the surface is given by $μ^{N_{C}} v_{i}^{N_{C}} = {\bar{v}}_{i}$ . Thus, it follows that Optimization B

\bar{D} (σ, N_{E}, N_{C}, c^{*}) := \min_{μ, R, L \geq 0} F (c^{*} ∥ \bar{v})

is equivalent to (6). Since $v$ is explicitly known as a function of $(R, L)$ , Optimization B (21) is a more tractable optimization problem than (6). Note that in this setting the function $D (σ, N_{E}, N_{C}, c^{*})$ in (21) is independent of the rates $μ$ . See Appendix 4—figure 1 for a flow chart of the two optimization schemes.

Appendix 4—figure 1

Download asset Open asset

Flow chart showing the optimization schemes for Optimization A and B.

We prove that $D_{\min}^{A} = D_{\min}^{B}$ by showing the set of all $c$ is equal to the set of all $v$ . We additionally establish that the optimum $v_{\min} = c_{\min}$ .

While 21 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 physicochemical 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 to define constrained Optimization B:

\bar{D} (σ, N_{C}, N_{E}, c^{*}) := \begin{aligned} min_{μ, R, L} & F (c^{*} ‖ \bar{v}) \\ s.t. & R_{min} \leq R (j, k, α) \leq R_{max}, \\ μ_{min} \leq μ^{(j)} \leq μ_{max} \\ 1 \leq ℓ_{α}^{(j)} \leq N_{S} . \end{aligned}

The upper and lower bounds on the rates $R$ and $μ$ are estimated in Appendix 2: $μ_{\max} = 1$ /min (resp. $μ_{\min} = .01$ /min) and $R_{\max} = 20$ /min (resp. $R_{\min} = .018$ /min).

Appendix 5

Numerical scheme for performing the non-convex optimization

We solve Optimization C using the numerical scheme detailed below. The optimization problem consists of minimizing a non-convex objective with linear box constraints. We use the MATLAB FMINCON function to solve this optimization. We use SQP, a gradient-based iterative optimization scheme for solving optimizations with non-linear differentiable objective and constraints. Since our problem is non-convex and SQP only gives local minima, we initialize the algorithm with many random initial points. We use SOBOLSET function of MATLAB to generate space filling pseudo-random numbers. We have taken 1000 initializations for each $N_{E}, N_{C}$ , and $σ$ value. We have taken 50 equally spaced points between 0 and 1 to explore the $σ$ -space for Figure 3. Some minor fluctuations in $D$ due to non-convexity of the objective function in the final results were smoothed out by taking the convex hull of the $D$ vs. $σ$ graph. The results for $σ_{m i n} (N_{E}, N_{C})$ and $D (σ_{m i n}, N_{E}, N_{C})$ (Figure 4) were obtained by adding $σ$ to the optimization vector and then performing the optimization.

A similar numerical scheme was used to optimize diversity.

Appendix 6

Analytical solution when $N_{s} ≫ 1$

It is possible to obtain analytical expressions for the steady-state glycan distribution in the limit $N_{s} ≫ 1$ when the glycan index $k$ can be approximated by a continuous variable. In this case, (15)–(16) can be cast as differential equations,

\begin{array}{ll} \frac{d c_{k}^{(1)}}{d k} \approx c_{k}^{(1)} - c_{k - 1}^{(1)} \\ = (\frac{\sum_{α = 1}^{N_{E}} R (1, k - 1, α) \exp (- σ | k - 1 - l_{α}^{(1)} |)}{μ^{(1)} + \sum_{α = 1}^{N_{E}} R (1, k, α) \exp (- σ | k - l_{α}^{(1)} |)} - 1) c_{k - 1}^{(1)} \\ \approx - (\frac{μ^{(1)} + \frac{d}{d k} \sum_{α = 1}^{N_{E}} R (1, k, α) \exp (- σ | k - l_{α}^{(1)} |)}{μ^{(1)} + \sum_{α = 1}^{N_{E}} R (1, k, α) \exp (- σ | k - l_{α}^{(1)} |)}) c_{k}^{(1)}, \end{array}

and

\begin{array}{ll} \frac{d c_{k}^{(j)}}{d k} \approx c_{k}^{(j)} - c_{k - 1}^{(j)} \\ = \frac{μ^{(j - 1)}}{μ^{(j)} + \sum_{α = 1}^{N_{E}} R (j, k, α) \exp (- σ | k - l_{α}^{(j)} |)} c_{k}^{(j - 1)} \\ - (\frac{μ^{(j)} + \frac{d}{d k} \sum_{α = 1}^{N_{E}} R (j, k, α) \exp (- σ | k - l_{α}^{(j)} |)}{μ^{(j)} + \sum_{α = 1}^{N_{E}} R (j, k, α) \exp (- σ | k - l_{α}^{(j)} |)}) c_{k}^{(j)} \end{array}

for $j = 2, \dots, N_{C}$ . In equation 22 and equation 23,

\begin{array}{ll} \frac{d}{d k} \sum_{α = 1}^{N_{E}} R (j, k, α) \exp (- σ | k - l_{α}^{(j)} |) \\ = \sum_{α = 1}^{N_{E}} R (j, k, α) σ \exp (- σ | k - l_{α}^{(j)} |) (1 - 2 I (k \geq l_{α})) + R^{'} (j, k, α) \exp (- σ | k - l_{α}^{(j)} |) \end{array}

where the indicator function $I (\cdot)$ is equal to 1 if the argument is true, and 0 otherwise and $R^{'} (j, k, α)$ is the derivative of $R (j, k, α)$ with respect to $k$ .

Define a vector function $C (k) \in ℝ_{c}^{N}$ of the continuous variable $k$ by $C (k) = [c_{k}^{(1)}, c_{k}^{(2)}, \dots c_{k}^{(N_{C})}]$ . Then, (equation 22) and (equation 23) can be written as

\frac{d C (k)}{d k} = M (k) C (k)

where the matrix $M (k)$ is given by

M (k) = [\begin{matrix} A^{(1)} (k) & 0 & 0 & 0 & \dots 0 \\ B^{(2)} (k) & A^{(2)} (k) & 0 & 0 & \dots 0 \\ 0 & B^{(3)} (k) & A^{(3)} (k) & 0 & \dots 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & \dots & 0 & B^{(N_{C})} (k) & A^{(N_{C})} (k) \end{matrix}]

with

\begin{array}{ll} A^{(j)} (k) = - \frac{μ^{(j)} + \frac{d}{d k} \sum_{α = 1}^{N_{E}} R (j, k, α) \exp (- σ | k - l_{α}^{(j)} |)}{μ^{(j)} + \sum_{α = 1}^{N_{E}} R (j, k, α) \exp (- σ | k - l_{α}^{(j)} |)} \\ B^{(j)} (k) = \frac{μ^{(j - 1)}}{μ^{(j)} + \sum_{α = 1}^{N_{E}} R (j, k, α) \exp (- σ | k - l_{α}^{(j)} |)} \end{array}

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 (25)

C (k) = C_{0} \exp (Ω (k))

is written in terms of the Magnus function $Ω (k) = \sum_{n = 1}^{\infty} Ω (n, k)$ , obtained from the Baker–Campbell–Hausdorff formula (Blanes et al., 2009),

\begin{array}{ll} Ω (1, k) & = \int_{0}^{k} M (k_{1}) d k_{1} \\ Ω (2, k) & = \frac{1}{2} \int_{0}^{k} d k_{1} \int_{0}^{k_{1}} d k_{2} [M (k_{1}), M (k_{2})] \\ Ω (3, k) & = \frac{1}{6} \int_{0}^{k} d k_{1} \int_{0}^{k_{1}} d k_{2} \int_{0}^{k_{2}} d k_{3} [M (k_{1}), [M (k_{2}), M (k_{3})]] + [M (k_{3}), [M (k_{2}), M (k_{1})]] \\ \dots \dots \end{array}

where $\dots$ is the commutator, and the higher order terms in … contain higher order nested commutators.

Here, we establish conditions under which the series $\sum_{n = 1}^{\infty} Ω (n, k)$ that defines solution $C (k)$ to the differential equation (25) converges. We also solve (25) for some special cases.

The commutator

[M (k_{1}), M (k_{2})] = [\begin{matrix} 0 & 0 & 0 & 0 & \dots & 0 \\ a_{21} & 0 & 0 & 0 & \dots & 0 \\ a_{31} & a_{32} & 0 & 0 & \dots & 0 \\ 0 & a_{42} & a_{43} & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & \dots & a_{n, n - 2} & a_{n, n - 1} & 0 \end{matrix}]

where

\begin{array}{ll} a_{i, i - 1} = A^{(i - 1)} (k_{2}) B^{(i)} (k_{1}) + A^{(i)} (k_{1}) B^{(i)} (k_{2}) - A^{(i - 1)} (k_{1}) B^{(i)} (k_{2}) + A^{(i)} (k_{2}) B^{(i)} (k_{1}) \\ a_{i, i - 2} = B^{(i - 1)} (k_{2}) B^{(i)} (k_{1}) - B^{(i - 1)} (k_{1}) B^{(i)} (k_{2}) \end{array}

The general form of $Ω (n, k)$ is given by Blanes et al., 2009

Ω (n, k) = \frac{z_{n}}{n!} \int_{0}^{k} d k_{1} \int_{0}^{k_{1}} d k_{2} \dots \int_{0}^{k_{n - 2}} d k_{n - 1} \int_{0}^{k_{n - 1}} d k_{n} \sum_{l} W_{l} M (k_{p_{1}^{l}}) M (k_{p_{2}^{l}}) \dots M (k_{p_{n}^{l}})

where $(p_{1}^{(l)}, p_{2}^{(l)} \dots p_{n}^{(l)})$ is a permutation of $(1, 2, 3, \dots n)$ , $W_{l} \in {- 1, 1}$ , and $z_{n} \in 1, \dots n$ .

Let $\bar{A} = \max_{k, l, m} | M_{l, m} (k) |$ . Define

\bar{M} = [\begin{matrix} \bar{A} & 0 & 0 & 0 \dots & 0 \\ \bar{A} & \bar{A} & 0 & 0 \dots & 0 \\ 0 & \bar{A} & \bar{A} & 0 \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & \dots & 0 & \bar{A} & \bar{A} \end{matrix}]

We can bound all the matrix elements of $Ω (n, k)$ in the following way:

\begin{array}{ll} Ω_{l m} (n, k) & \leq z_{n} {\bar{M}}_{l, m}^{n} \int_{0}^{k} d k_{1} \int_{0}^{k_{1}} d k_{2} \dots \int_{0}^{k_{n - 1}} d k_{n} \\ = z_{n} {\bar{M}}^{n} |_{l m} \frac{k^{n}}{n!} \end{array}

The matrix

{\bar{M}}^{n} = [\begin{matrix} a_{11} & 0 & 0 & 0 & \dots & 0 \\ a_{21} & a_{22} & 0 & 0 & \dots & 0 \\ a_{31} & a_{32} & a_{33} & 0 & \dots & 0 \\ a_{41} & a_{42} & a_{43} & a_{44} & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ a_{n 1} & \dots & a_{n, n - 2} & a_{n, n - 1} & a_{n n} \end{matrix}]

where $a_{l m} = S_{l m} (n) {\bar{A}}^{n}$ for appropriately defined polynomials $S_{l, m} (n)$ . Thus, it follows that $Ω_{l m} \leq z_{n} S_{l m} (n) {(A^{*})}^{n} \frac{k^{n}}{n!}$ and $Ω_{l, m} (k) \leq \sum_{n = 1}^{\infty} z_{n} S_{l, m} (n) {(A^{*})}^{n} \frac{k^{n}}{n!}$ . Consequently, the series will converge if $\bar{A} k < 1$ , that is, $k \leq \frac{1}{\bar{A}}$ . Assuming $μ^{(j)} = μ \forall j$ , we can bound $\bar{A}$ as

\bar{A} \leq \max_{j, k} (\frac{μ + σ \sum_{α = 1}^{N_{E}} R (j, k, α) \exp (- σ | k - l_{α}^{(j)} |)}{μ + \sum_{α = 1}^{N_{E}} R (j, k, α) \exp (- σ | k - l_{α}^{(j)} |)} + \frac{\sum_{α = 1}^{N_{E}} R^{'} (j, k, α) \exp (- σ | k - l_{α}^{(j)} |)}{μ + \sum_{α = 1}^{N_{E}} R (j, k, α) \exp (- σ | k - l_{α}^{(j)} |)})

Since the parameters $μ, σ, R (j, k, α), l_{α}^{(j)}$ , and $N_{E}$ are finite and positive, and $R^{'} (j, k, α)$ is finite, $\bar{A}$ has a finite upper bound, implying that $k$ is always greater than zero, and the series has a 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: (i) $(N_{E} = 1, N_{C} = 1)$ and (ii) $(N_{E} = 1, N_{C} = 2)$ .

i. $N_{E} = 1, N_{C} = 1$ : The solution of the differential equation is given by

c (k) = {\begin{cases} c_{0} e^{- k} {(\frac{μ + R \exp (- σ (l - k))}{μ + R \exp (- σ l)})}^{(1 / σ) - 1} & k \leq l \\ c (l) e^{- (k - l)} {(\frac{μ + R}{μ + R \exp (- σ (k - l))})}^{(1 / σ) + 1} & k > l \end{cases}

A representative concentration profile is plotted in Appendix 6—figure 1. 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 6—figure 1

Download asset Open asset

Glycan concentration profile calculated from the model using (a) formula (31) for $N_{E} = N_{C} = 1$ and (b) formulae (32)–(36) for $N_{E} = 1, N_{C} = 2$ .

ii. $N_{E} = 1, N_{C} = 2$ : The concentration profile $c^{(2)}$ in cisterna-2 can be obtained from the following calculation. Let $l^{(j)}$ denote the ‘length’ of the enzyme in cisterna $j = 1, 2$ . For $k \leq \min {l^{(1)}, l^{(2)}}$ ,

\begin{array}{ll} c^{(2)} (k) = & c_{0} μ^{(1)} e^{- k} (\frac{μ^{(2)} + R^{(2)} \exp (- σ (l^{(2)} - k))}{μ^{(1)} + R^{(1)} e^{- σ l^{(1)}}})^{(1 / σ) - 1} \int_{0}^{k} \frac{(μ^{(1)} + R^{(1)} \exp (- σ (l^{(1)} - k)))^{(1 / σ) - 1}}{(μ^{(2)} + R^{(2)} \exp (- σ (l^{(2)} - k)))^{1 / σ}} d k \\ + c^{(2)} (0) e^{- k} {(\frac{μ^{(2)} + R^{(2)} e^{- σ (l^{(2)} - k)}}{μ^{(2)} + R^{(2)} e^{- σ l^{(2)}}})}^{(1 / σ) - 1} \end{array}

Next, consider the case where $l^{(1)} \leq l^{(2)}$ . Then, for $l^{(1)} < k \leq l^{(2)}$

\begin{array}{ll} c^{(2)} (k) & = c^{(1)} (l^{(1)}) μ^{(1)} e^{- (k - l^{(1)})} (μ^{(1)} + R^{(1)})^{(1 / σ) + 1} (μ^{(2)} + R^{(2)} \exp (- σ (l^{(2)} - k)))^{(1 / σ) - 1} \\ \int_{l^{(1)}}^{k} \frac{(μ^{(2)} + R^{(2)} \exp (- σ (l^{(2)} - k)))^{- 1 / σ}}{(μ^{(1)} + R^{(1)} \exp (- σ (k - l^{(1)})))^{(1 / σ) + 1}} d k \\ + c^{(2)} (l^{(1)}) e^{- (k - l^{(1)})} {(\frac{μ^{(2)} + R^{(2)} e^{- σ (l^{(2)} - k)}}{μ^{(2)} + R^{(2)} e^{- σ (l^{(2)} - l^{(1)})}})}^{(1 / σ) - 1} \end{array}

and for $l^{(1)} \leq l^{(2)} < k$ ,

\begin{array}{ll} c^{(2)} (k) & = c^{(1)} (l^{(1)}) μ^{(1)} e^{- (k - l^{(1)})} {(\frac{μ^{(1)} + R^{(1)}}{μ^{(2)} + R^{(2)} \exp (- σ (k - l^{(2)}))})}^{(1 / σ) + 1} \\ \int_{l^{(2)}}^{k} \frac{(μ^{(2)} + R^{(2)} \exp (- σ (k - l^{(2)})))^{1 / σ}}{(μ^{(1)} + R^{(1)} \exp (- σ (k - l^{(1)})))^{(1 / σ) + 1}} d k \\ + c^{(2)} (l^{(2)}) e^{- (k - l^{(2)})} {(\frac{μ^{(2)} + R^{(2)}}{μ^{(2)} + R^{(2)} e^{- σ (k - l^{(2)})}})}^{(1 / σ) + 1} \end{array}

Next, the case where $l^{(1)} \geq l^{(2)}$ . For $l^{(2)} < k \leq l^{(1)}$ ,

\begin{array}{ll} c^{(2)} (k) & = c_{0} μ^{(1)} e^{- k} \frac{(μ^{(1)} + R^{(1)} e^{- σ l^{(1)}})^{1 - (1 / σ)}}{(μ^{(2)} + R^{(2)} \exp (- σ (k - l^{(2)})))^{(1 / σ) + 1}} \int_{l^{(2)}}^{k} \frac{(μ^{(1)} + R^{(1)} \exp (- σ (l^{(1)} - k)))^{(1 / σ) - 1}}{(μ^{(2)} + R^{(2)} \exp (- σ (k - l^{(2)})))^{- 1 / σ}} d k \\ + c^{(2)} (l^{(2)}) e^{l^{(2)} - k} {(\frac{μ^{(2)} + R^{(2)}}{μ^{(2)} + R^{(2)} e^{- σ (k - l^{(2)})}})}^{(1 / σ) + 1} \end{array}

For $^{(2)} \leq l^{(1)} < k$ ,

\begin{array}{ll} c^{(2)} (k) & = c^{(1)} (l^{(1)}) μ^{(1)} e^{- (k - l^{(1)})} {(\frac{μ^{(1)} + R^{(1)}}{μ^{(2)} + R^{(2)} \exp (- σ (k - l^{(2)}))})}^{(1 / σ) + 1} \\ \int_{l^{(2)}}^{k} \frac{(μ^{(2)} + R^{(2)} \exp (- σ (k - l^{(2)})))^{1 / σ}}{(μ^{(1)} + R^{(1)} \exp (- σ (k - l^{(1)})))^{(1 / σ) + 1}} d k \\ + c^{(2)} (l^{(1)}) e^{- (k - l^{(1)})} {(\frac{μ^{(2)} + R^{(2)} e^{- σ (l^{(1)} - l^{(2)})}}{μ^{(2)} + R^{(2)} e^{- σ (k - l^{(2)})}})}^{(1 / σ) + 1} \end{array}

The integrals in (32)–(36) can be evaluated numerically. The result of the numerical computation is shown in Appendix 6—figure 1.

Appendix 6—figure 2a–d plots the glycan profile $c_{k} v s . 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}$ (Appendix 6—figure 2a and 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}$ (Appendix 6—figure 2b–d).
Keeping everything else fixed, decreasing enzyme specificity increases the spread of the distribution around the peaks (Appendix 6—figure 2a and c).

Appendix 6—figure 2

Download asset Open asset

Glycan profile ${c_{k} : k = 1, \dots, N_{s}}$ as a function of specificity $σ$ (a, c) and reaction rates $R$ (b, d).

(a) $N_{E} = N_{C} = 1, (R = 50, μ = 1, l = 10)$ . c_k decreases exponentially with $k$ for very low and very high $σ$ ; however, the decay rate is lower at low $σ$ . For intermediate values of $σ$ , the distribution has *exactly* two peaks, one of which is at $k = 0$ , and eventually decays exponentially. The width of the distribution is a decreasing function of $σ$ . (b) $N_{E} = N_{C} = 1, (σ = 0.1, μ = 1, l = 10)$ . At low $R$ , c_k is concentrated at low $k$ . The proportion of higher index glycans in an increasing function of $R$ . (c) $N_{E} = N_{C} = 2, (R = 40, μ = 1, [l_{1}^{(1)}, l_{2}^{(1)}, l_{1}^{(2)}, l_{2}^{(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. (d) $N_{E} = N_{C} = 2, (R = 40, μ = 1, [l_{1}^{(1)}, l_{2}^{(1)}, l_{1}^{(2)}, l_{2}^{(2)}] = [10, 30, 50, 70])$ . As in the plots in (b), increasing $R$ shifts the peaks towards higher index glycans and the proportion of higher index glycan increases.

Appendix 7

Extended distortion model shows lack of apparent symmetry

The results in Figure 4 seem to imply that there is an approximate $N_{E} - N_{C}$ symmetry in the model, that is, increasing either $N_{E}$ or $N_{C}$ affects the fidelity, optimal enzyme specificity, and the sensitivity in approximately the same way. This is a consequence of the distortion model we are using for calculating the binding probabilities of substrates with enzymes, which allows every enzyme $α$ to in principle catalyze any reaction in any cisterna. This allowed for the ideal enzyme length $l_{α}^{(j)}$ in Equation 3 to vary across the cisternae in an unconstrained manner, leading to simplification in the calculation. We now find that by changing this aspect of the model the apparent symmetry between $N_{E} - N_{C}$ is lifted. A more reasonable model for the ideal enzyme length is given by

ℓ_{α}^{(j)} = ℓ_{α}^{(0)} + δ ℓ_{α}^{(j)}, δ ℓ_{α}^{(j)} \in [- ℓ_{b}^{(j)}, ℓ_{b}^{(j)}],

that is, the nominal length $ℓ_{α}^{(0)}$ can be distorted in a cisterna by a correction $δ ℓ_{α}^{(j)}$ but within a specified bound $ℓ_{b}^{(j)}$ that is not subject to optimization. One can render some enzymes inactive in certain cisternae by choosing appropriate values of $ℓ_{α}^{(0)}$ and $δ ℓ_{α}^{(j)}$ . For small values for the bound $ℓ_{b}$ , for example, $l_{b} / N_{s} \leq 0.2$ (here, $N_{s} - 1$ is the number of enzymatic reactions), the decrease in $\bar{D}$ on increasing $N_{C}$ is small compared to increasing $N_{E}$ (see Appendix 7—figure 1). On the other hand, for large $ℓ_{b}$ , for example, $l_{b} / N_{S} \geq 0.3$ , there is an approximate symmetry between $N_{E}$ and $N_{C}$ (see Appendix 7—figure 1). Here, we have taken the bounds to be compartment independent, that is, $ℓ_{b}^{(j)} = ℓ_{b}$ .

Appendix 7—figure 1

Download asset Open asset

Optimum fidelity ${\bar{D}}_{K L}$ as a function of $(N_{E}, N_{C})$ for different values of $ℓ_{b} / N_{s}$ , where $ℓ_{b}$ bounds the deformation in the ideal length $ℓ_{α}^{(0)}$ of an enzyme $α = 1, \dots, N_{E}$ .

Small values of $ℓ_{b}$ restrict all enzymes from working in all cisternae and all substrates, where large value of $ℓ_{b}$ removes this constraint.

Appendix 8

Redundancies and non-convexity of the optimization

Validation of the numerical optimization scheme

In order to test whether our numerical optimization procedure is able to converge to the global minimum, we run the following test. We generate 100 random values of $(μ, R, L, σ)$ within their respective ranges for a problem instance with $(N_{E} = 2, N_{C} = 2)$ . The sampled value for $(μ, R, L, σ)$ is used to generate concentration profiles that are then used as the target distribution for the optimization. Since the target distribution is achievable, the optimal value of the constrained Optimization B for these sampled targets is $\bar{D} = 0$ . We solve the constrained Optimization B using our numerical scheme. The average optimal value $\bar{D}$ across all sampled values was 9.1835e-07, 30 out of 100 values were exactly zero, and the highest $\bar{D}$ was 1.1761e-05. Therefore, the optimization scheme was able to recover the concentration profiles almost exactly. Next, we ask whether the optimization problem recovers the value of $(μ, R, L, σ)$ that was used to create the particular target distribution. We were able to recover $σ$ exactly, except in cases where the concentration profile was almost a delta function at the first glycan (see Appendix 8—figure 1). This is because $σ$ decides the typical width of the empirical distribution, and hence the optimal $σ$ is determined by the typical width of the target distribution, except in the pathological case of a concentration profile that is almost a delta function at the first glycan – such a concentration profile can be made produced for any value $σ$ by simply making transport $μ$ very fast as compared to the reaction rates.

Appendix 8—figure 1

Download asset Open asset

Recovering the $σ$ values for different target distribution.

Note that barring four data points all other optimized $σ$ values (red dots) exactly overlap with the corresponding target $σ$ (diamonds).

We note that the optimization in $(μ, R, L)$ is not convex and leads to many equally good minimas corresponding to different values of $(μ, R, L)$ . The resulting redundancies in the model and their importance are discussed next.

Degeneracy in the model

Recall that in equation 7, we defined

R_{e f f} (j, k) = \frac{\sum_{α} R_{α}^{(j)} \exp (- σ | k - l_{α}^{(j)} |)}{μ^{(j)}}

In terms of these renormalized rates, the steady-state glycan concentration can be written as

{\bar{c}}_{k}^{(j)} = \frac{R_{e f f} (j, k) {\bar{c}}_{k - 1}^{(j)} + {\bar{c}}_{k}^{j - 1}}{1 + R_{e f f} (j, k)},

that is, the concentration is only a function of $R_{e f f} (k, j)$ . Thus, any combination of $(μ, R, L, σ)$ that maps to the same value of $R_{e f f}$ will result in the same concentration profile and will be indistinguishable from the perspective of the objective function. Additionally, the mapping from $R_{e f f}$ to the concentration profile $\bar{c}$ also has degeneracy. We show these redundancies in the schematic below, which shows a systematic reduction in dimension to 1 (scalar), which is the quantity we optimize,

\begin{array}{ccccccc} R, L, μ, σ & ⟶ & R_{e f f} & ⟶ & \bar{c} & ⟶ & F \\ N_{C} + 2 N_{E} N_{C} + 1 & N_{C} (N_{S} - 1) & N_{S} - 1 & 1 \end{array}

Since $F (c_{T} | | \bar{c}) = 0$ if and only if $c_{T} = \bar{c}$ , it follows that the last mapping does not have redundancy. Some of the sources of degeneracies in the mapping from $(μ, R, L, σ)$ to $R_{e f f}$ are as follows:

For fixed $(σ, L)$ , setting $R_{α}^{(j)} \leftarrow γ R_{α}^{(j)}$ and $μ_{j} \leftarrow γ^{- 1} μ_{j}$ leaves $R_{e f f}$ invariant.
Permutations in the $α$ index leave $R_{e f f}$ invariant. Thus, there are at least ${(N_{E}!)}^{N_{C}}$ distinct minima that map to the same value of $R_{e f f}$ , and therefore, the same concentration $\bar{c}$ .

Additionally, there are degeneracies coming from the optimization which depend on the target distribution $c_{T}$ . Having discussed the sources of degeneracies of the optimized solution, we now discuss the distribution of the optimized solutions.

Distribution of minima

To study the behaviour of the optimization algorithm for different initial points, we numerically investigate the distribution of function values at different local minima. Since the dimension of the optimization problem is $N_{C} + 2 N_{E} N_{C} + 1$ , which is large, we divide the optimization space into a grid of $I = n_{p}^{(N_{C} + 2 N_{E} N_{C} + 1)}$ points. We did this numerical experiment for $(N_{E}, N_{C}) \in {(1, 1), (1, 2), (2, 1), (2, 2)}$ . The value of $n_{p} = 3$ for $(N_{E}, N_{C}) = (2, 2)$ and $n_{p} = 4$ for the rest. The target distribution for all the cases is a single Gaussian with mean 20, standard deviation 5, with support on $1 \leq k \leq 20$ . The results of this numerical experiment are summarized in Appendix 8—figure 2 and Appendix 8—table 1, from which we deduce the following:

A large fraction of the initial starting points converge to a set of degenerate minima with objective function value exactly equal to the global minimum. These minima are a result of the degeneracies of the optimization problem.
There are other local minima with objective value very close to (but not equal) the global minimum. Most initial points converge to one of these two sets of minima.
Finally, there are a small set of local minima with significantly higher objective values. These correspond to minima with $σ = 0$ . The fraction of initial points that converge to such minima reduces as the dimension of the optimization space increases.

Appendix 8—table 1

Distribution of local minima.

$N_{E}$	$N_{C}$	$\min {\bar{D}}_{K L}$	$\max {\bar{D}}_{K L}$	Fraction of initial conditions within ${\bar{D}}_{K L} \leq 0.0228$
$N_{E}$	$N_{C}$	$\min {\bar{D}}_{K L}$	$\max {\bar{D}}_{K L}$		1	1	0.0228	0.44	0.56
2	1	0.0081	0.44	0.73
1	2	0.0051	0.29	0.70
2	2	1.17e-4	0.29	0.84

Appendix 8—figure 2

Download asset Open asset

$\bar{D}$ for various initial conditions, sorted in increasing order for clarity.

This clearly shows the fraction of initial conditions for which the optimized $\bar{D}$ is small (see Appendix 8—table 1).

Appendix 9

Robustness of optimization to small perturbations

We analyse the sensitivity to small perturbations around the optimal point by calculating the Hessian of the KL divergence,

H (i, j) = {\frac{\partial^{2}}{\partial X_{i} \partial X_{j}} F |}_{X_{min}}

Here, $X = [μ, R, L, σ]$ denotes the entire set of optimization variables (note that the entries in $X$ are normalized by their respective range and do not carry physical dimensions). We calculated the eigenvalues, denoted by $λ_{i}$ , and eigenvectors, denoted by $V_{i}$ , of the Hessian matrix to identify the stiff and sloppy directions (Gutenkunst et al., 2007; Machta et al., 2013) in the optimization space. The eigenvectors of the Hessian matrix can be grouped in $R, L, μ$ , and $σ$ directions by looking for the most dominant component in the eigenvector. We find that most of the eigenvectors have significant entries along the direction of only one of the optimization variables $μ, R, L, σ$ , for example, in Figure 6a, the eigenvectors 21–36 have significant entries only in the $L$ directions. There is, however, a small number of eigenvectors that have entries over more than one optimization direction, for example, the eigenvector with $σ$ dominant direction has some $μ$ component as well (Figure 6a).

Stiff and sloppy directions

We find that the eigenvalues of the eigenvectors dominated by $σ$ and some $μ$ , $L$ directions are orders of magnitude higher than for those dominated by the $R$ directions (see Figure 6b). This suggests that $\bar{D}$ has a valley-like structure around the optimal, with $R$ and some $L$ being the flat or sloppy directions.

The fact that enzyme specificity $σ$ and some of the $L$ directions are stiff should not be surprising since the typical width and position of peaks in the synthesized distribution are primarily controlled by $σ$ and $L$ . We have already shown that $\bar{D}$ is a sharp convex function of $σ$ for low values of $(N_{E}, N_{C})$ (see Figure 3), which gradually flattens out as we increase $(N_{E}, N_{C})$ .

The fact that transport rate $μ$ is a stiff direction is surprising! The stiffness in $μ$ is due to the fact that the optimal $μ$ is always at the lower bound, and with even slight increase in $μ$ , the transport becomes too fast for the reactions to be able to produce the intermediate products. For the $(R, L)$ -dominated eigenvectors, there are bands of sloppy direction and stiff directions. We define the average stiffness in $μ, R, L$ , and $σ$ by a weighted average of eigenvalues, where the weight is given by the strength of the corresponding components of the eigenvector.

⟨ λ ⟩_{μ} = \ln (\sum_{i} w_{i}^{(μ)} λ_{i}), ⟨ λ ⟩_{R} = \ln (\sum_{i} w_{i}^{(R)} λ_{i}), ⟨ λ ⟩_{L} = \ln (\sum_{i} w_{i}^{(L)} λ_{i}), ⟨ λ ⟩_{σ} = \ln (\sum_{i} w_{i}^{(σ)} λ_{i})

Here, $w_{i}^{(μ)} = \sum_{j \in μ} | V_{i, j} | / \sum_{j} | V_{i, j} |$ , $w_{i}^{(R)} = \sum_{j \in R} | V_{i, j} | / \sum_{j} | V_{i, j} |$ , $w_{i}^{(L)} = \sum_{j \in L} | V_{i, j} | / \sum_{j} | V_{i, j} |$ and $w_{i}^{(μ)} = \sum_{j \in σ} | V_{i, j} | / \sum_{j} | V_{i, j} |$ .

Figure 6c shows ${⟨ λ ⟩}_{μ}$ , ${⟨ λ ⟩}_{R}$ , ${⟨ λ ⟩}_{L}$ , and ${⟨ λ ⟩}_{σ}$ as a function of $N_{C}$ for fixed $N_{E} = 4$ . The average stiffness in $R$ directions, ${⟨ λ ⟩}_{R}$ , is considerably lower than the average stiffness in $σ$ , $μ$ and $L$ directions. $σ$ is the stiffest direction but the stiffness decreases on increasing the $N_{C}$ . Interestingly, the stiffness along $L$ directions increases on increasing $N_{C}$ .

We now define the total average stiffness $⟨ λ ⟩ = \log (\frac{\sum λ_{i}}{N_{C} + 2 N_{E} N_{C} + 1})$ , that is, log of the sum of the eigenvalues divided by the dimension of the optimization problem, in the space of $N_{E}, N_{C}$ . We find that the average stiffness is higher for low values of $(N_{E}, N_{C})$ as compared to higher values of $(N_{E}, N_{C})$ , with a few exceptions; and eventually, the average stiffness settles to a fixed low value (Figure 6d).

Appendix 10

Insensitivity of diversity to threshold c_th

Since the threshold $c_{t h}$ used to count the number of glycan species in the definition of diversity is arbitrary (see section ‘Strategies to achieve high glycan diversity’), we here show that the qualitative results we obtain are independent of this choice (Appendix 10—figure 1).

Appendix 10—figure 1

Download asset Open asset

Diversity vs. $N_{C}$ for different values of $σ$ keeping $N_{E} = 1$ fixed, for three different values of the threshold, $c_{t h} = \frac{1}{N_{s}}, \frac{1}{2 N_{s}}, \frac{1}{4 N_{s}}$ .

Changing the value of the threshold $c_{t h}$ only changes the saturation value of the diversity curve.

Appendix 11

Table of symbols used in the article (Appendix 11—table 1)

Appendix 11—table 1

Table of symbols and their definitions.

Symbol	Definition
$c_{k}^{(j)}$	Concentration of $k$ th glycan in $j$ th compartment
$μ^{(j)}$	Transport rate from $j$ th to $j + 1$ compartment
$σ$	Specificity of the enzymes
$ℓ_{α}^{(j)}$	Ideal substrate length for αth enzyme in $j$ th compartment
$M (j, k, α)$	Enzyme parameter related to the Michaelis constant $K_{M}$
$V (j, k, α)$	Enzyme parameter related to $V_{m a x}$
$R (j, α)$	Reaction parameter for Optimization B
$D (c^{*} ∥ c)$	KL divergence between $c^{*}$ and $c$
$F (c^{*} \| \| c)$	Fidelity, KL divergence normalized by the entropy of the target $c^{*}$
$\bar{D} (σ, N_{E}, N_{C}, c^{*})$	$m i n_{μ, R, L} F (c^{*} \| \| c)$
$R_{e f f} (j, k)$	Effective reaction rate of $k$ th reaction in $j$ th compartment

Data availability

The current manuscript is a computational study, so no data have been generated for this manuscript. The following repository on github contains the code and the data (numerical data + Mass Spec data) that are used in the paper: https://github.com/alkeshyadav/Glycosylation, (copy archived at URL swh:1:rev:a46c6eb76c5f07458d07e44267f48bbaaff6fc5a).

The following data sets were generated

1. Yadav A
(2022) GitHub
ID Glycosylation. Glycan processing in the Golgi: optimal information coding and constraints on cisternal number and enzyme specificity.

https://github.com/alkeshyadav/Glycosylation

References

Book
1. Alberts B
(2002)
Molecular Biology of the Cell

Garland Science.
- Google Scholar
(2013) Evolving views of pectin biosynthesis
Annual Review of Plant Biology 64:747–779.
https://doi.org/10.1146/annurev-arplant-042811-105534
- PubMed
- Google Scholar
1. Bacharoglou A
(2010) Approximation of probability distributions by convex mixtures of Gaussian measures
Proceedings of the American Mathematical Society 138:2619–2628.
https://doi.org/10.1090/S0002-9939-10-10340-2
- Google Scholar
1. Bar-Even A
2. Milo R
3. Noor E
4. Tawfik DS
(2015) The Moderately Efficient Enzyme: Futile Encounters and Enzyme Floppiness
Biochemistry 54:4969–4977.
https://doi.org/10.1021/acs.biochem.5b00621
- PubMed
- Google Scholar
1. Bard F
2. Chia J
(2016) Cracking the Glycome Encoder: Signaling, Trafficking, and Glycosylation
Trends in Cell Biology 26:379–388.
https://doi.org/10.1016/j.tcb.2015.12.004
- PubMed
- Google Scholar
1. Becker B
2. Melkonian M
(1996) The secretory pathway of protists: spatial and functional organization and evolution
Microbiological Reviews 60:697–721.
https://doi.org/10.1128/mr.60.4.697-721.1996
- Google Scholar
1. Berninsone PM
2. Hirschberg CB
(2000) Nucleotide sugar transporters of the Golgi apparatus
Current Opinion in Structural Biology 10:542–547.
https://doi.org/10.1016/s0959-440x(00)00128-7
- PubMed
- Google Scholar
1. Blanes S
2. Casas F
3. Oteo JA
4. Ros J
(2009) The Magnus expansion and some of its applications
Physics Reports 470:151–238.
https://doi.org/10.1016/j.physrep.2008.11.001
- Google Scholar
1. Bonner JT
(1998) The origins of multicellularity
Integrative Biology 1:27–36.
https://doi.org/10.1002/(SICI)1520-6602(1998)1:1<27::AID-INBI4>3.0.CO;2-6
- Google Scholar
Book
1. Boyd S
2. Vandenberghe L
(2004) Convex Optimization
Cambridge University Press.
https://doi.org/10.1017/CBO9780511804441
- Google Scholar
1. Burton RA
2. Fincher GB
(2014) Evolution and development of cell walls in cereal grains
Frontiers in Plant Science 5:456.
https://doi.org/10.3389/fpls.2014.00456
- PubMed
- Google Scholar
1. Caffaro CE
2. Hirschberg CB
(2006) Nucleotide sugar transporters of the Golgi apparatus: from basic science to diseases
Accounts of Chemical Research 39:805–812.
https://doi.org/10.1021/ar0400239
- PubMed
- Google Scholar
1. Carroll SB
(2001) Chance and necessity: the evolution of morphological complexity and diversity
Nature 409:1102–1109.
https://doi.org/10.1038/35059227
- PubMed
- Google Scholar
(2010) Sensors and regulators of intracellular pH
Nature Reviews. Molecular Cell Biology 11:50–61.
https://doi.org/10.1038/nrm2820
- PubMed
- Google Scholar
1. Changeux JP
2. Edelstein SJ
(2005) Allosteric mechanisms of signal transduction
Science (New York, N.Y.) 308:1424–1428.
https://doi.org/10.1126/science.1108595
- PubMed
- Google Scholar
Book
1. Cover TM
2. Thomas JA
(2012)
Elements of Information Theory

John Wiley & Sons.
- Google Scholar
1. Cummings RD
2. Pierce JM
(2014) The challenge and promise of glycomics
Chemistry & Biology 21:1–15.
https://doi.org/10.1016/j.chembiol.2013.12.010
- PubMed
- Google Scholar
Data
1. Cummings RD
2. Crocker P
(authors) (2020) Functional Glycomics Database
Consortium for Functional Glycomics.

http://www.functionalglycomics.org
(2001) Negative regulation of T-cell activation and autoimmunity by Mgat5 N-glycosylation
Nature 409:733–739.
https://doi.org/10.1038/35055582
- PubMed
- Google Scholar
(2009) Metabolism, cell surface organization, and disease
Cell 139:1229–1241.
https://doi.org/10.1016/j.cell.2009.12.008
- PubMed
- Google Scholar
1. Dmitrieff S
2. Rao M
3. Sens P
(2013) Quantitative analysis of intra-Golgi transport shows intercisternal exchange for all cargo
PNAS 110:15692–15697.
https://doi.org/10.1073/pnas.1303358110
- PubMed
- Google Scholar
(2007) Identification and characterization of COPIa- and COPIb-type vesicle classes associated with plant and algal Golgi
PNAS 104:163–168.
https://doi.org/10.1073/pnas.0609818104
- PubMed
- Google Scholar
1. Drickamer K
2. Taylor ME
(1998) Evolving views of protein glycosylation
Trends in Biochemical Sciences 23:321–324.
https://doi.org/10.1016/s0968-0004(98)01246-8
- PubMed
- Google Scholar
1. Du P
2. Stolovitzky G
3. Horvatovich P
4. Bischoff R
5. Lim J
6. Suits F
(2008) A noise model for mass spectrometry based proteomics
Bioinformatics (Oxford, England) 24:1070–1077.
https://doi.org/10.1093/bioinformatics/btn078
- PubMed
- Google Scholar
1. Dwek RA
(1996) Glycobiology: Toward Understanding the Function of Sugars
Chemical Reviews 96:683–720.
https://doi.org/10.1021/cr940283b
- PubMed
- Google Scholar
(2013) Glycosphingolipids: synthesis and functions
The FEBS Journal 280:6338–6353.
https://doi.org/10.1111/febs.12559
- PubMed
- Google Scholar
(2018) Cell type transcriptome atlas for the planarian Schmidtea mediterranea
Science (New York, N.Y.) 360:eaaq1736.
https://doi.org/10.1126/science.aaq1736
- PubMed
- Google Scholar
1. Fisher P
2. Ungar D
(2016) Bridging the Gap between Glycosylation and Vesicle Traffic
Frontiers in Cell and Developmental Biology 4:15.
https://doi.org/10.3389/fcell.2016.00015
- PubMed
- Google Scholar
1. Fisher P
2. Spencer H
3. Thomas-Oates J
4. Wood AJ
5. Ungar D
(2019) Modeling Glycan Processing Reveals Golgi-Enzyme Homeostasis upon Trafficking Defects and Cellular Differentiation
Cell Reports 27:1231–1243.
https://doi.org/10.1016/j.celrep.2019.03.107
- PubMed
- Google Scholar
1. Free SJ
(2013) Fungal cell wall organization and biosynthesis
Advances in Genetics 81:33–82.
https://doi.org/10.1016/B978-0-12-407677-8.00002-6
- PubMed
- Google Scholar
1. Gabius HJ
(2018) The sugar code: Why glycans are so important
Bio Systems 164:102–111.
https://doi.org/10.1016/j.biosystems.2017.07.003
- PubMed
- Google Scholar
1. Gagneux P
2. Varki A
(1999) Evolutionary considerations in relating oligosaccharide diversity to biological function
Glycobiology 9:747–755.
https://doi.org/10.1093/glycob/9.8.747
- PubMed
- Google Scholar
1. Glick BS
2. Malhotra V
(1998) The curious status of the Golgi apparatus
Cell 95:883–889.
https://doi.org/10.1016/s0092-8674(00)81713-4
- PubMed
- Google Scholar
1. Gow NAR
2. Hube B
(2012) Importance of the Candida albicans cell wall during commensalism and infection
Current Opinion in Microbiology 15:406–412.
https://doi.org/10.1016/j.mib.2012.04.005
- PubMed
- Google Scholar
(2007) Universally sloppy parameter sensitivities in systems biology models
PLOS Computational Biology 3:0030189.
https://doi.org/10.1371/journal.pcbi.0030189
- PubMed
- Google Scholar
1. Han X
2. Zhou Z
3. Fei L
4. Sun H
5. Wang R
6. Chen Y
7. Chen H
8. Wang J
9. Tang H
10. Ge W
11. Zhou Y
12. Ye F
13. Jiang M
14. Wu J
15. Xiao Y
16. Jia X
17. Zhang T
18. Ma X
19. Zhang Q
20. Bai X
21. Lai S
22. Yu C
23. Zhu L
24. Lin R
25. Gao Y
26. Wang M
27. Wu Y
28. Zhang J
29. Zhan R
30. Zhu S
31. Hu H
32. Wang C
33. Chen M
34. Huang H
35. Liang T
36. Chen J
37. Wang W
38. Zhang D
39. Guo G
(2020) Construction of a human cell landscape at single-cell level
Nature 581:303–309.
https://doi.org/10.1038/s41586-020-2157-4
- PubMed
- Google Scholar
1. Hayashi T
2. Kaida R
(2011) Functions of xyloglucan in plant cells
Molecular Plant 4:17–24.
https://doi.org/10.1093/mp/ssq063
- PubMed
- Google Scholar
(1998) Transporters of nucleotide sugars, ATP, and nucleotide sulfate in the endoplasmic reticulum and Golgi apparatus
Annual Review of Biochemistry 67:49–69.
https://doi.org/10.1146/annurev.biochem.67.1.49
- PubMed
- Google Scholar
(2007) Systems analysis of N-glycan processing in mammalian cells
PLOS ONE 2:e713.
https://doi.org/10.1371/journal.pone.0000713
- PubMed
- Google Scholar
Preprint
1. Jaiman A
2. Thattai M
(2018) Algorithmic Biosynthesis of Eukaryotic Glycans
bioRxiv.
https://doi.org/10.1101/440792
- Google Scholar
1. Kellokumpu S
(2019) Golgi pH, Ion and Redox Homeostasis: How Much Do They Really Matter?
Frontiers in Cell and Developmental Biology 7:93.
https://doi.org/10.3389/fcell.2019.00093
- PubMed
- Google Scholar
Book
1. Kirschner MW
2. Gerhart JC
(2008)
The Plausibility of Life

Yale University Press.
- Google Scholar
(2015) Keeping track of the growing number of biological functions of chitin and its interaction partners in biomedical research
Glycobiology 25:469–482.
https://doi.org/10.1093/glycob/cwv005
- PubMed
- Google Scholar
1. Krambeck FJ
2. Betenbaugh MJ
(2005) A mathematical model of N-linked glycosylation
Biotechnology and Bioengineering 92:711–728.
https://doi.org/10.1002/bit.20645
- PubMed
- Google Scholar
1. Krambeck FJ
2. Bennun SV
3. Narang S
4. Choi S
5. Yarema KJ
6. Betenbaugh MJ
(2009) A mathematical model to derive N-glycan structures and cellular enzyme activities from mass spectrometric data
Glycobiology 19:1163–1175.
https://doi.org/10.1093/glycob/cwp081
- PubMed
- Google Scholar
1. Kumar P
2. Yang M
3. Haynes BC
4. Skowyra ML
5. Doering TL
(2011) Emerging themes in cryptococcal capsule synthesis
Current Opinion in Structural Biology 21:597–602.
https://doi.org/10.1016/j.sbi.2011.08.006
- PubMed
- Google Scholar
(2002) Structure of the Golgi and distribution of reporter molecules at 20 degrees C reveals the complexity of the exit compartments
Molecular Biology of the Cell 13:2810–2825.
https://doi.org/10.1091/mbc.01-12-0593
- PubMed
- Google Scholar
(1998) Measurement of cytosolic, mitochondrial, and Golgi pH in single living cells with green fluorescent proteins
PNAS 95:6803–6808.
https://doi.org/10.1073/pnas.95.12.6803
- PubMed
- Google Scholar
(2013) Parameter space compression underlies emergent theories and predictive models
Science (New York, N.Y.) 342:604–607.
https://doi.org/10.1126/science.1238723
- PubMed
- Google Scholar
Book
1. MacKay DJ
(2003)
Information Theory, Inference and Learning Algorithms

Cambridge University Press.
- Google Scholar
(2014) The Cell Biology of Cellulose Synthesis
Annual Review of Plant Biology 65:69–94.
https://doi.org/10.1146/annurev-arplant-050213-040240
- PubMed
- Google Scholar
(2017) Different Golgi ultrastructure across species and tissues: Implications under functional and pathological conditions, and an attempt at classification
Tissue and Cell 49:186–201.
https://doi.org/10.1016/j.tice.2016.12.002
- PubMed
- Google Scholar
(2003) Tomographic evidence for continuous turnover of Golgi cisternae in Pichia pastoris
Molecular Biology of the Cell 14:2277–2291.
https://doi.org/10.1091/mbc.e02-10-0697
- PubMed
- Google Scholar
(1965) On the nature of allosteric transitions: a plausible model
Journal of Molecular Biology 12:88–118.
https://doi.org/10.1016/s0022-2836(65)80285-6
- PubMed
- Google Scholar
1. Moremen KW
2. Haltiwanger RS
(2019) Emerging structural insights into glycosyltransferase-mediated synthesis of glycans
Nature Chemical Biology 15:853–864.
https://doi.org/10.1038/s41589-019-0350-2
- PubMed
- Google Scholar
1. Nam H
2. Lewis NE
3. Lerman JA
4. Lee DH
5. Chang RL
6. Kim D
7. Palsson BO
(2012) Network context and selection in the evolution to enzyme specificity
Science (New York, N.Y.) 337:1101–1104.
https://doi.org/10.1126/science.1216861
- PubMed
- Google Scholar
(2004) Rhamnogalacturonan II: structure and function of a borate cross-linked cell wall pectic polysaccharide
Annual Review of Plant Biology 55:109–139.
https://doi.org/10.1146/annurev.arplant.55.031903.141750
- PubMed
- Google Scholar
1. Pauly M
2. Gille S
3. Liu L
4. Mansoori N
5. de Souza A
6. Schultink A
7. Xiong G
(2013) Hemicellulose biosynthesis
Planta 238:627–642.
https://doi.org/10.1007/s00425-013-1921-1
- PubMed
- Google Scholar
1. Peracchi A
(2018) The Limits of Enzyme Specificity and the Evolution of Metabolism
Trends in Biochemical Sciences 43:984–996.
https://doi.org/10.1016/j.tibs.2018.09.015
- PubMed
- Google Scholar
(2019) Translation of genome to glycome: role of the Golgi apparatus
FEBS Letters 593:2390–2411.
https://doi.org/10.1002/1873-3468.13541
- PubMed
- Google Scholar
Book
1. Price N.
2. Stevens L.
(1999)
Fundamentals of Enzymology: The cell and molecular biology of catalytic proteins

Oxford University Press.
- Google Scholar
1. Roseman S
(2001) Reflections on glycobiology
The Journal of Biological Chemistry 276:41527–41542.
https://doi.org/10.1074/jbc.R100053200
- PubMed
- Google Scholar
1. Sachdeva H
2. Barma M
3. Rao M
(2016) Nonequilibrium description of de novo biogenesis and transport through Golgi-like cisternae
Scientific Reports 6:1–10.
https://doi.org/10.1038/srep38840
- PubMed
- Google Scholar
(2014) Unique, polyfucosylated glycan-receptor interactions are essential for regeneration of Hydra magnipapillata
ACS Chemical Biology 9:147–155.
https://doi.org/10.1021/cb400486t
- PubMed
- Google Scholar
1. Savir Y
2. Tlusty T
(2007) Conformational proofreading: the impact of conformational changes on the specificity of molecular recognition
PLOS ONE 2:e468.
https://doi.org/10.1371/journal.pone.0000468
- PubMed
- Google Scholar
1. Sengupta D
2. Linstedt AD
(2011) Control of organelle size: the Golgi complex
Annual Review of Cell and Developmental Biology 27:57–77.
https://doi.org/10.1146/annurev-cellbio-100109-104003
- PubMed
- Google Scholar
1. Sens P
2. Rao M
(2013) (Re)modeling the Golgi
Methods in Cell Biology 118:299–310.
https://doi.org/10.1016/B978-0-12-417164-0.00018-5
- PubMed
- Google Scholar
(2019) Stem cell differentiation trajectories in Hydra resolved at single-cell resolution
Science (New York, N.Y.) 365:eaav9314.
https://doi.org/10.1126/science.aav9314
- PubMed
- Google Scholar
(2018) Identification of multiple isomeric core chitobiose-modified high-mannose and paucimannose N-glycans in the planarian Schmidtea mediterranea
The Journal of Biological Chemistry 293:6707–6720.
https://doi.org/10.1074/jbc.RA117.000782
- PubMed
- Google Scholar
1. Umaña P
2. Bailey JE
(1997) A mathematical model of N-linked glycoform biosynthesis
Biotechnology and Bioengineering 55:890–908.
https://doi.org/10.1002/(SICI)1097-0290(19970920)55:6<890::AID-BIT7>3.0.CO;2-B
- PubMed
- Google Scholar
(1983) Terminal alpha (1 leads to 4)-linked N-acetylglucosamine: A characteristic constituent of duodenal-gland mucous glycoproteins in rat and pig. A high-resolution 1H-NMR study
Biochimica et Biophysica Acta 747:107–116.
https://doi.org/10.1016/0167-4838(83)90128-0
- PubMed
- Google Scholar
1. Varki A
(1998) Factors controlling the glycosylation potential of the Golgi apparatus
Trends in Cell Biology 8:34–40.
https://doi.org/10.1016/s0962-8924(97)01198-7
- PubMed
- Google Scholar
Book
1. Varki A
(2009)
Essentials of Glycobiology

Cold Spring Harbor Laboratory Press.
- Google Scholar
1. Varki A
(2011) Evolutionary forces shaping the Golgi glycosylation machinery: why cell surface glycans are universal to living cells
Cold Spring Harbor Perspectives in Biology 3:a005462.
https://doi.org/10.1101/cshperspect.a005462
- PubMed
- Google Scholar
1. Varki A
(2017) Biological roles of glycans
Glycobiology 27:3–49.
https://doi.org/10.1093/glycob/cww086
- PubMed
- Google Scholar
1. Wills C
2. Green DR
(1995) A genetic herd-immunity model for the maintenance of MHC polymorphism
Immunological Reviews 143:263–292.
https://doi.org/10.1111/j.1600-065x.1995.tb00679.x
- PubMed
- Google Scholar
1. Winterburn PJ
2. Phelps CF
(1972) The significance of glycosylated proteins
Nature 236:147–151.
https://doi.org/10.1038/236147a0
- PubMed
- Google Scholar
1. Yang M
2. Fehl C
3. Lees KV
4. Lim EK
5. Offen WA
6. Davies GJ
7. Bowles DJ
8. Davidson MG
9. Roberts SJ
10. Davis BG
(2018) Functional and informatics analysis enables glycosyltransferase activity prediction
Nature Chemical Biology 14:1109–1117.
https://doi.org/10.1038/s41589-018-0154-9
- PubMed
- Google Scholar

Decision letter

Felix Campelo

Reviewing Editor; The Barcelona Institute of Science and Technology, Spain
Suzanne R Pfeffer

Senior Editor; Stanford University School of Medicine, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Decision letter after peer review:

[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]

Thank you for submitting your work entitled "Glycan processing in the Golgi – optimal information coding and constraints on cisternal number and enzyme specificity" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and a Senior Editor. The reviewers have opted to remain anonymous.

This manuscript presents a clever approach to study from a theoretical point of view how glycosylation reactions occurring at the Golgi membranes define a cell-specific glycan profile. In particular, by optimizing a metric, the authors find that larger number of enzymes and/or compartments (Golgi cisternae) are required for increasingly complex glycan profiles. In addition, they found that enzymes should not be very specific nor too sloppy. This manuscript touches upon an important topic in cell biology (origin of glycan diversity and link to Golgi complex architecture), which has been relatively unexplored in the past. However, although there are no major concerns about the validity of the model presented in the paper (see separate reviews for more detailed comments), there is a general consensus among the reviewers that there are some key points that this manuscript should resolve before opting for publication at eLife. Since these points will most likely require an extensive revision of the manuscript, we decided to proceed with a rejection at this moment. We encourage you, if you wish, to resubmit a stronger new paper in the future; or to opt for submission elsewhere for its timely publication.

The main concerns raised by the reviewers are the following:

1) It is not clear what new biology we learn from the results of this paper (see e.g. Rev. 1's report). It's been suggested that the reason for stacking is to give proteins enough time to undergo proper glycosylation and that cisternae number correlates with the complexity of glycans. Since the paper does not include any new experimental evidence on how the complexity of cellular glycans leads to more cisternae in the Golgi stack, it would be very valuable for a biological audience if the authors could present a list of model predictions that could be experimentally tested.

2) For the results of the model to be more robust and probably to also learn a bit more about cell specific glycobiology, a more systematic study of available glycan profiles would very much improve the story (see in particular, Rev. 2, point 2; and also Rev. 3, point 3).

3) The manuscript is very complex for non theorists, since it is presented as a rather technical manuscript. It is well written and organized for specialists. However, it is not clear whether this manuscript, in the current form, will be of interest for a broad biological audience. Along these lines, the conclusions and discussion might need to be toned down, given the large number of simplifications of the model (which is fine for a theoretical audience, but could be over-assumed by a non-specialist biological audience) (see Rev. 3, point 1).

Reviewer #1:

The manuscript by Yadav et al. uses a theoretical model of glycan processing to investigate the constraints placed on the enzymatic modification process in order to achieve the complex final glycan distribution. The authors introduce a multi-compartment, many-species model of chemical modifications and develop an optimization algorithm to minimize the divergence between the resulting glycan distribution and that observed from data on immune cells. The authors find that a minimum number of enzymes and reactions chambers (cisternae) are necessary to recapitulate the complex features of the experimental distribution. They also find that enzymes should have an optimal specificity: they should not be too specific to only one particular substrate, nor should they be too "sloppy" by interacting with all substrates with similar rates. Finally, the authors elucidate control parameters of their model that maximize the diversity of glycans.

I am torn about this study because while the approach is elegant and principled and the results sensible and interpretable, the scope is rather technical and therefore I am skeptical as to whether this will have the broad audience and impact that is expected for eLife. The work combines elements from several fields (cell biology, chemical modeling, computational optimization, information theory), which is commendable, well-presented, and appropriate for an interdisciplinary audience. At the same time, however, I found myself wondering what we really learn at the end of the day. That the components must be complex (sufficient number of enzymes and compartments) in order for the result to be complex (many-peaked glycan distribution) seems unavoidable. That the specificity must be moderate also seems expected in order to not be limited by a one-to-one matching between enzyme and glycan product, but also not produce a broad, sloppy distribution. Part of the problem may also be that it is unclear what contribution is made in terms of predictions for new experiments (which is often a necessary criterion for a high-impact theoretical study): apart from using actual data as an optimization target, the study largely recapitulates and rationalizes existing features of the glycan machinery, rather than suggesting perturbations for, or placing quantitative limits on, new experimental investigations. My inclination is therefore to suspect that this work may be more suitable for a more specialized journal, e.g. on cell biology or chemical physics.

Reviewer #2:

The manuscript "Glycan processing in the Golgi – optimal information coding and constraints on cisternal number and enzyme specificity" by Yadav et al. presents a model for glycan production by the Golgi cisternae. Glycans serve various cellular functions, such as markers for cell identity, and it is important that a function-specific profile (distribution) of glycans is obtained at the end of the production process. The authors show how to choose the reaction constants in their model to reach a target profile (distribution) of glycans. They discuss how the adequacy between the target and model profiles depend on the physico-chemical features of the reactants.

While the problem studied by the authors is interesting and their analysis reports valuable qualitative results, e.g. on the consequence of enzyme specificity on the profile statistics, I list below various points that I think would require further studies and discussions:

1. The writing of the manuscript is quite technical in some sections, in particular around equation (5) and (6), and involves superfluous notations: what is the added value in introducing V, M, L with respect to the already defined omega, P, …? I would recommend the authors to simplify the writing and really emphasize what is important here. Several definitions or assumptions remain vague and not discussed. For instance, equation (3) define the binding probability between enzyme α and substrate k in terms of "shape" vectors. On line 250 these shape vectors are defined as scalars, in particular the shape of enzyme k is identified with the index itself. Why is this a good choice (apart from reducing the dimension of the parameter space), and how could one test the validity this apparently quite strong hypothesis?

2. The manuscript is mostly theoretical with limited connections with experimental data and validations. It would be nice to have a systematic study of available glycan profiles and comparison of the corresponding inferred model parameters. In addition, it is necessary to validate the inference/optimization procedure on synthetic profiles generated by stochastic simulations of some prescribed network, to check that the correct network parameters are correctly recovered (or some other solutions giving back the same profile, see point 4 below).

3. While the KL divergence (7) is a reasonable measure of dissimilarity between two distributions the relevance of the ratio (9) on page 10 is not at all clear to me. More precisely, I do not see why the entropy of the target distribution is the relevant quantity here. Rather the KL divergence between the model and target distribution could be compared to the mean value of the KL divergence between any two glycan profiles resulting from the noisy production steps.

4. Since the function of the minimize over, see eq. (8) is non convex, there may exist many local minima. The authors report in Appendix 8 they perform grid-like search. How close are those minima in terms of the objective function? This question is important as the model studied here seems to be overparametrized, in particular when N_c is large. It may also be relevant in terms of evolution. More precisely, how were some glycan profiles and the corresponding biochemical networks selected in the course of evolution? On top of selection pressure to "orthogonalize" the profiles is it possible that some distributions are (approximately) realized by more networks than others and were therefore entropically favored?

Reviewer #3:

The paper by Yadav et al. contributes to an important and largely unexplored topic in Life Sciences, namely, to the understanding of glycosylation. The authors introduce a mathematical model of glycosylation in the Golgi apparatus and use the model to investigate how the complexity (diversity) and fidelity of the plasma membrane glycan distribution depends on parameters such as the number of Golgi cisternae or enzyme specificity. The paper is well written and makes the effort to present a rather complex topic in an accessible way by leaving some of the details in the appendices.

After defining the quantitative measures for complexity and fidelity of glycan distributions, the paper builds a general model of glycosylation in the Golgi by assuming Michaelis-Menten type reactions and by taking into account parameters such as the number, specificity and distribution of enzymes across cisternae. However, the parameter space becomes so large that, before solving it numerically, the authors need to reduce it by making a number of simplifications (they linearize the distortion energy; assign one enzyme specificity σ to all reactions …). The authors then ask the question how different parameters influence the outcome and what are the optimal parameters for a given target glycan distribution. Finally, the authors explore the strategies to achieve high glycan diversity.

My comments are listed below:

1. The authors should be more careful (modest?) with the generalization of their results and conclusions. The results are obtained by a very simplified model and it is not obvious that they can be generalized.

For example, in the Abstract they state: "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." However, this is only true within the simplified model. Maybe, complex distributions could be obtained even with a small number of cisternae if one allowed biochemistry that is more complex. The statement above would be more accurate if it was turned around, e.g., "we show that multiple cisternae and precise enzyme partitioning in the Golgi can lead to complex glycan distributions."

The readers of a mathematical journal would very well understand that the results are not general but should be interpreted within the model proposed. On the other hand, the readers of eLife could be misled by statements which are too general. This is even more true for this paper, which starts with a rather general and complex model and then simplifies it on the go. While the authors do summarize the simplifications in the Discussion, an inexperienced reader could still miss them and take the results too literally. I would suggest that the authors amend the whole paper according to this comment to avoid possible misinterpretations.

2. An interesting result of the modeling is that the role of parameters NE and NC can be rather symmetrical – the fidelity of glycan distribution can be improved in a similar way by increasing either of these parameters (Figure 4). Why is that? Is this a consequence of a symmetrical role of these parameters in the reduced (simplified) parameter space, or can this symmetry be generalized and can have a biological merit?

3. The paper mentions examples of cells that have a complex Golgi and a complex glycan distribution. Is the opposite also true, e.g., that cells with a simple Golgi do not have complex glycan distributions? Or can a complex Golgi sometimes result in a low glycan complexity?

4. One wonders how robust the glycosylation process is according to the model proposed? Figure 5b neatly shows that large enzyme repartitioning can lead to a completely different glycan distribution. But how do small variations of the optimal reaction distribution affect the outcome? Would it be possible to quantify the robustness and analyze it systematically (not just show a few random examples)? How is robustness affected by the number of cisternae and enzymes? What are biological implications of this result?

5. Figure 5b shows that a complex Golgi can also lead to a simple glycan distribution. This implies that having a large NC and NE is not sufficient for a complex glycan distribution. Could the authors discuss this? What are additional requirements that are needed for a complex glycan distribution?

6. What is "likelihood," which is used in some of the graphs? Is it related to the minimum achievable Kullback-Leibler (KL) divergence according to Kullback-Leibler metric defined in Eq. 7 (and used in Figures 3, 4 and 5)? If so, the authors could use only one measure of likelihood/divergence in all the graphs.

7. Why does the mouse glycan distribution have a lower "likelihood" than the human one, even with an increasing number of GMM (Figure 1b)? Is this a coincidence or can this be generalized?

8. The paper would be easier to read if all the parameters were listed in one place, and if one letter was not used to denote different things (e.g., "k" denotes complexity, but it is also used as the index for glycosylation reactions).

9. Is the final distribution of glycans on PM presented on the right in Figure 2 conceptually the same as the glycan distributions presented in Figure 1? If so, use the same notation (and axes labels) everywhere.

https://doi.org/10.7554/eLife.76757.sa1

Author response

[Editors’ note: the authors resubmitted a revised version of the paper for consideration. What follows is the authors’ response to the first round of review.]

Reviewer #1:

The manuscript by Yadav et al. uses a theoretical model of glycan processing to investigate the constraints placed on the enzymatic modification process in order to achieve the complex final glycan distribution. The authors introduce a multi-compartment, many-species model of chemical modifications and develop an optimization algorithm to minimize the divergence between the resulting glycan distribution and that observed from data on immune cells. The authors find that a minimum number of enzymes and reactions chambers (cisternae) are necessary to recapitulate the complex features of the experimental distribution. They also find that enzymes should have an optimal specificity: they should not be too specific to only one particular substrate, nor should they be too "sloppy" by interacting with all substrates with similar rates. Finally, the authors elucidate control parameters of their model that maximize the diversity of glycans.

I am torn about this study because while the approach is elegant and principled and the results sensible and interpretable, the scope is rather technical and therefore I am skeptical as to whether this will have the broad audience and impact that is expected for eLife. The work combines elements from several fields (cell biology, chemical modeling, computational optimization, information theory), which is commendable, well-presented, and appropriate for an interdisciplinary audience. At the same time, however, I found myself wondering what we really learn at the end of the day. That the components must be complex (sufficient number of enzymes and compartments) in order for the result to be complex (many-peaked glycan distribution) seems unavoidable. That the specificity must be moderate also seems expected in order to not be limited by a one-to-one matching between enzyme and glycan product, but also not produce a broad, sloppy distribution. Part of the problem may also be that it is unclear what contribution is made in terms of predictions for new experiments (which is often a necessary criterion for a high-impact theoretical study): apart from using actual data as an optimization target, the study largely recapitulates and rationalizes existing features of the glycan machinery, rather than suggesting perturbations for, or placing quantitative limits on, new experimental investigations. My inclination is therefore to suspect that this work may be more suitable for a more specialized journal, e.g. on cell biology or chemical physics.

We are happy that the referee finds the paper “elegant”, “commendable” and “appropriate for an interdisciplinary audience”. Despite this the referee does not feel that the manuscript provides new predictions or suggestions for new experiments and simply “rationalizes existing features of the glycan machinery”.

While it is true that our study provides rationalisations for existing (but not well appreciated) features of the cellular glycan machinery, it does so from a fresh perspective, namely from a distinct information theoretic point of view. This should not be considered as a problem, since one of the roles of theory is to view disparate observations through the lens of a general framework. Our abstract information theoretic considerations are a novel conceptual way of approaching the biochemical synthesis of information carrying molecules, and make quantitative, the purely verbal descriptions in, 2. We hope such a perspective would prove useful in the future.

In addition, our theoretical study also makes testable predictions, provides suggestions for future experiments and opens up a scope for new enquiry. Our rewriting of the manuscript emphasises this aspect – we do this by (i) including and analysing more data from different species, that together highlight the relation between organismal/glycan complexity and the complexity of the Golgi machinery; and (ii) highlighting the predictions and suggestions that follow from this study (this has involved some more calculations). While the abstraction and generality of the description helps in bringing together several aspects and highlighting new questions, it can at best make qualitative predictions. Our task therefore has been to to get at a qualitative understanding using quantitative methods and thereby to arrive at general principles. This approach has been highly successful in several branches of physics, a rather famous example is the Landau theory of phase transitions, and has been passionately, and we think correctly, advocated for Biology in [3, 4].

We now highlight the new conceptual contributions, predictions and offshoots/implications of our theoretical study and suggestions for future experiments -

1. Since an important function of the glycan spectrum is cell type/niche identification, it seems natural to relate glycan complexity to organismal complexity taken to be associated with the number of cell types in the organism [5, 6]. Here, we provide a measure of the complexity of the glycan distribution of a given cell type using MSMS data. Using this we have now analysed the MSMS data from Hydra, Planaria, Mammalian cells.

2. We have constructed a generic and calculable model of chemical synthesis of the glycan spectrum (molecular information) consistent with the phenomenology of enzymatic reactions and cellular compartmentalization. This represents a calculable model, as opposed to a black-box model, of the Golgi machinery.

3. Constructing a high fidelity representation of a complex target distribution, such as those observed in real cells, requires a complex Golgi machinery with multiple cisternae, precise enzyme partitioning and control on enzyme specificity. This definition of fidelity of the glycan code, allows us to provide a quantitative argument for the evolutionary requirement of multiple-compartments. While it is possible to produce complex glycan distributions in one compartment using a large number of enzymes, such a design would inevitably require a more elaborate genetic cost.

4. One of the predictions of the model is that we should see a qualitative trend of increasing Golgi compartments on increasing the glycan complexity, keeping everything else fixed. This can be empirically verified by looking at the glycan profiles, ultrastructure of Golgi and the number of glycosylation enzymes across many species.

5. The model solution is degenerate, in the sense that there many equally good global minima. These degeneracies are both continuous and discrete. The continuous degeneracies corresponds to regions in the reaction rate – transport rate space, moving along which will not change the concentration profile, thus ensuring robustness to internal noise. This suggests that the distribution is robust to slight cell-to-cell variations in these kinetic parameters.

6. Our model implies that close to a local minima the inter-cisternal transport rate µ and the specificity of the enzymes σ are stiff directions, i.e. the cell should exercise tighter control on µ and σ as compared to the reaction rates. The reaction rates close to the local minima are sloppy directions, moving along these directions does not change the glycan profile much.

7. 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. 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. Sloppy or promiscuous enzymes bring in the potential for evolvability [7], sloppiness allows the system to be stable to random mutations in proteins or variations in the target distribution.

8. Taken together, our quantitative analysis of the tradeoffs has profound implications for nonequilibrium self assembly of the Golgi cisternae, suggesting that the nonequilibrium control of cisternal number must involve a coupling of nonequilibrium self assembly of cisternae with enzymatic chemical reaction kinetics [8].

Reviewer #2:

The manuscript "Glycan processing in the Golgi – optimal information coding and constraints on cisternal number and enzyme specificity" by Yadav et al. presents a model for glycan production by the Golgi cisternae. Glycans serve various cellular functions, such as markers for cell identity, and it is important that a function-specific profile (distribution) of glycans is obtained at the end of the production process. The authors show how to choose the reaction constants in their model to reach a target profile (distribution) of glycans. They discuss how the adequacy between the target and model profiles depend on the physico-chemical features of the reactants.

While the problem studied by the authors is interesting and their analysis reports valuable qualitative results, e.g. on the consequence of enzyme specificity on the profile statistics, I list below various points that I think would require further studies and discussions:

1. The writing of the manuscript is quite technical in some sections, in particular around equation (5) and (6), and involves superfluous notations: what is the added value in introducing V, M, L with respect to the already defined omega, P, …? I would recommend the authors to simplify the writing and really emphasize what is important here.

We thank the referee for this suggestion. V and M are physically measurable macroscopic quantities (vectors) related to the Michaelis constant and maximum velocity V_max of the enzymes, for which we place bounds that we take from literature (see Appendix 2). These are defined in terms of the microscopic rates ω_f,ω_b and ω_c, see equation 4 of the manuscript. However to make the manuscript clearer we have now relegated the expression for the concentrations in terms of M,V and L (equations 5 and 6 of the manuscript) to Appendix 3. We have also simplified the writing as recommended, highlighting essential features.

Several definitions or assumptions remain vague and not discussed. For instance, equation (3) defines the binding probability between enzyme α and substrate k in terms of "shape" vectors. On line 250 these shape vectors are defined as scalars, in particular the shape of enzyme k is identified with the index itself. Why is this a good choice (apart from reducing the dimension of the parameter space), and how could one test the validity this apparently quite strong hypothesis?

We now state all the assumptions that we are making at the beginning of the Results section. In general, the shape of an object (enzyme) can be described by a shape tensor akin to a moment of inertia tensor (or a fabric tensor that expands the shape in spherical harmonics). Restricting the shape deformations to a scalar is tantamount to a single mode approximation. This is essentially the approximation used when a single reaction coordinate is taken to represent a multidimensional reaction landscape.

In a chemical reaction between macromolecules, the reaction centres are often buried in the molecule and not easily accessible, unless reactants have the appropriate shape (lock-key) or can be induced to have the appropriate shape (induced fit). Thus, reactions rates depend on the shape and its deformability, which may be modelled by a general shape tensor (we call this distortion). A simple model for distortion is a scalar model, as used for instance in [9]. We have chosen to parametrize the optimal length of the enzyme by a integer multiple of a basic unit. This choice merely simplifies the analysis. Normally in the analysis of chemical kinetics, this would be subsumed in the rates, since the experimentally measured rates are a function of the shape distortion, pH etc. Since in our work, we are interested in analysing the tradeoffs involving enzymatic parameters, we need to provide a “microscopic” model of enzyme reactions that allows us to proceed with our optimization analysis in a tractable manner.

2. The manuscript is mostly theoretical with limited connections with experimental data and validations. It would be nice to have a systematic study of available glycan profiles and comparison of the corresponding inferred model parameters.

What the referee is asking for is precisely the goal of glycan engineering (see Discussion in the revised manuscript). This is a mature and serious field of study, and has contributed a lot to understanding the chemical profile of glycans in health and disease. There are a number of sequential enzymatic reaction models [10, 11, 12] that predict the N-glycan distribution based on the activities and levels of processing enzymes distributed in the Golgi-cisternae of mammalian cells, and make comparisons with N-glycan mass spectrum data. Models such as the KB2005 model [10, 11, 12] are extremely elaborate (with 22800 chemical reactions) and require very detailed chemical input parameters, which are either determined experimentally or obtained by fitting to the mass spectrometry data. These models have an important practical role to play, that of being able to predict the impact of the various chemical parameters on glycan distribution.

In contrast, here we pose a different question: Given an evolutionary drive for a specific cell type within a specific tissue microenvironment to maintain a target glycan distribution of fixed complexity, what are the tradeoffs in the chemical and physical parameters of glycan biosynthesis that realise this target.

This engenders a different strategy of analysis as advocated here, and allows us to get at a qualitative understanding using quantitative methods, as elaborated in the Response to Referee 1.

In addition, it is necessary to validate the inference/optimization procedure on synthetic profiles generated by stochastic simulations of some prescribed network, to check that the correct network parameters are correctly recovered (or some other solutions giving back the same profile, see point 4 below).

We agree with the referee and thank the referee for this suggestion. We have now added a new section in the manuscript (Please see Appendix 8 of the revised manuscript).

3. While the KL divergence (7) is a reasonable measure of dissimilarity between two distributions the relevance of the ratio (9) on page 10 is not at all clear to me. More precisely, I do not see why the entropy of the target distribution is the relevant quantity here. Rather the KL divergence between the model and target distribution could be compared to the mean value of the KL divergence between any two glycan profiles resulting from the noisy production steps.

The entropy of target distribution is a measure of the total statistical information in the distribution. Normalizing with respect to the entropy allows us to measure the information gained about the target distribution by approximating it with the empirical distribution, as the fraction of the total statistical information in the target distribution. This is a reasonable measure for comparing the effectiveness of the model for different target distributions. We have now made this point more explicit in the revised manuscript (Please see Equation 5 in section "Optimization Problem").

4. Since the function of the minimize over, see eq. (8) is non convex, there may exist many local minima. The authors report in Appendix 8 they perform grid-like search. How close are those minima in terms of the objective function? This question is important as the model studied here seems to be overparametrized, in particular when N_c is large. It may also be relevant in terms of evolution. More precisely, how were some glycan profiles and the corresponding biochemical networks selected in the course of evolution? On top of selection pressure to "orthogonalize" the profiles is it possible that some distributions are (approximately) realized by more networks than others and were therefore entropically favored?

The issue of parameter degeneracy and possible evolutionary implications is very interesting, and we thank the referee for bringing this up. We have conducted a detailed study, which we now report in Appendix 8 subsection “Degeneracy in the model”.

Reviewer #3:

The paper by Yadav et al. contributes to an important and largely unexplored topic in Life Sciences, namely, to the understanding of glycosylation. The authors introduce a mathematical model of glycosylation in the Golgi apparatus and use the model to investigate how the complexity (diversity) and fidelity of the plasma membrane glycan distribution depends on parameters such as the number of Golgi cisternae or enzyme specificity. The paper is well written and makes the effort to present a rather complex topic in an accessible way by leaving some of the details in the appendices.

After defining the quantitative measures for complexity and fidelity of glycan distributions, the paper builds a general model of glycosylation in the Golgi by assuming Michaelis-Menten type reactions and by taking into account parameters such as the number, specificity and distribution of enzymes across cisternae. However, the parameter space becomes so large that, before solving it numerically, the authors need to reduce it by making a number of simplifications (they linearize the distortion energy; assign one enzyme specificity σ to all reactions …). The authors then ask the question how different parameters influence the outcome and what are the optimal parameters for a given target glycan distribution. Finally, the authors explore the strategies to achieve high glycan diversity.

My comments are listed below:

1. The authors should be more careful (modest?) with the generalization of their results and conclusions. The results are obtained by a very simplified model and it is not obvious that they can be generalized.

For example, in the Abstract they state: "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." However, this is only true within the simplified model. Maybe, complex distributions could be obtained even with a small number of cisternae if one allowed biochemistry that is more complex. The statement above would be more accurate if it was turned around, e.g., "we show that multiple cisternae and precise enzyme partitioning in the Golgi can lead to complex glycan distributions."

The readers of a mathematical journal would very well understand that the results are not general but should be interpreted within the model proposed. On the other hand, the readers of eLife could be misled by statements which are too general. This is even more true for this paper, which starts with a rather general and complex model and then simplifies it on the go. While the authors do summarize the simplifications in the Discussion, an inexperienced reader could still miss them and take the results too literally. I would suggest that the authors amend the whole paper according to this comment to avoid possible misinterpretations.

We apologise for not being careful in stating our results and presenting the necessary caveats. We believe we have now worked on it to a considerable extent in the revised submission.

We find that even within the framework of our simple model it is possible to generate complex distributions in a small number of compartments by having a complex biochemistry (large number of processing enzymes). However, this complex biochemistry invokes an elaborate genetic cost on the cell. Therefore, within our framework the number of compartments required are both a function of the complexity of glycan distribution to be synthesized and the number of enzymes. This is a crucial insight and we have edited the manuscript to make this point and the underlying assumptions explicitly.

2. An interesting result of the modelling is that the role of parameters NE and NC can be rather symmetrical – the fidelity of glycan distribution can be improved in a similar way by increasing either of these parameters (Figure 4). Why is that? Is this a consequence of a symmetrical role of these parameters in the reduced (simplified) parameter space, or can this symmetry be generalized and can have a biological merit?

We thank the referee for bringing this up. This is an interesting observation. However, we should hasten to add that this apparent symmetry is a consequence of our simplified chemical and transport model, as we have verified by making numerous extensions of the model. We now discuss these issues in detail in the section "Tradeoffs between number of enzymes, number of cisternae and enzyme specificity to achieve given complexity" of the revised manuscipt.

3. The paper mentions examples of cells that have a complex Golgi and a complex glycan distribution. Is the opposite also true, e.g., that cells with a simple Golgi do not have complex glycan distributions? Or can a complex Golgi sometimes result in a low glycan complexity?

Since the complexity of the glycan distribution depends on both the complexity of the Golgi and the complexity of the chemical processing, this is a very difficult question to answer. We do feel that a systematic analysis of the enzyme chemistry, Golgi morphology and transport would shed more light on this issue, and indeed this is one of the suggestions that emerge from our analysis. For this paper, we have taken the glycan MS data of of a few cell types – planaria, hydra, mouse and human blood cells.

4. One wonders how robust the glycosylation process is according to the model proposed? Figure 5b neatly shows that large enzyme repartitioning can lead to a completely different glycan distribution. But how do small variations of the optimal reaction distribution affect the outcome? Would it be possible to quantify the robustness and analyze it systematically (not just show a few random examples)? How is robustness affected by the number of cisternae and enzymes? What are biological implications of this result?

Again, we thank the referee for bringing up this very interesting issue. We now provide a detailed analysis of the sensitivity to small perturbations around the optimal solution by calculating the Hessian of the KL divergence in the new section "Geometry of the Fidelity landscape" of the revised manuscript.

7. Why does the mouse glycan distribution have a lower "likelihood" than the human one, even with an increasing number of GMM (Figure 1b)? Is this a coincidence or can this be generalized?

Apart from updating the GMM procedure, we have also included data from two different organisms, Hydra and Planaria, in our study. Comparisons across these species tell a story when viewed from our framework.

Figure 5 compares the complexity of glycan profiles of Hydra, Planaria and Humans. Number of cell types in Hydra, Planaria and Humans are around 41 [15], 44 [16] and 103 [17] respectively, based on transcriptome analysis (these are lower bounds based on the main cell types, and especially for Planaria and Hydra, are subject to constant revision). These studies cluster thousands of cells from various tissues in the organism on the basis of statistical similarity of mRNA expression in the cell, with each cluster representing a cell type. This method is quite good for estimating the cell types in simpler organisms or single tissues but undercounts the number of cell types for humans due to limited cell number of each cell type in the sample. Based on morphology the number of cell types in humans is estimated to be around 200 [18].

Our analysis of the MSMS data of these organisms suggest that organism with fewer cell types have less complex glycan distribution.

It would be interesting to determine systematic correlations between complexity of Golgi, number of glycosylation enzymes and the complexity of the glycan distribution, in organisms such as algae (which typically possess many cisternae).

We have included this figure and a related discussion in Appendix 1 of the revised manuscript.

[1] Gabius, Hans-Joachim. ”The sugar code: Why glycans are so important.” Biosystems 164 (2018): 102111.

[2] Bard, Frederic, and Joanne Chia. ”Cracking the glycome encoder: signaling, trafficking, and glycosylation.” Trends in cell biology 26.5 (2016): 379-388.

[3] Sengupta, Anirvan M. Modeling biomolecular networks: an introduction to systems biology. Oxford University Press, 2008.

[4] W. Bialek, W., ”Biophysics: Searching for Principles”, (Princeton University Press, Princeton, 2002).

[5] Carroll, Sean B. ”Chance and necessity: the evolution of morphological complexity and diversity.” Nature 409.6823 (2001): 1102-1109.

[6] Bonner, John Tyler. ”The origins of multicellularity.” Integrative Biology: Issues, News, and Reviews: Published in Association with The Society for Integrative and Comparative Biology 1.1 (1998): 27-36.

[7] Kirschner, M.W. and Gerhart, J.C., ”The plausibility of life”, (Yale University Press, New Haven, 2005).

[8] Glick, B.S. and Malhotra, V., ”The curious status of the Golgi apparatus”, Cell 95: 883-889 (1998).

[9] Savir Y, Tlusty T. "Conformational proofreading: the impact of conformational changes on the specificity of molecular recognition". PloS one. 2007; 2(5):e468

[10] Umana, Pablo, and James E. Bailey. "A mathematical model of N linked glycoform biosynthesis." Biotechnology and bioengineering 55.6 (1997): 890-908.

[11] Krambeck, Frederick J., and Michael J. Betenbaugh. "A mathematical model of N linked glycosylation." Biotechnology and Bioengineering 92.6 (2005) 711-728.

[12] Krambeck, Frederick J., et al. "A mathematical model to derive N-glycan structures and cellular enzyme activities from mass spectrometric data." Glycobiology 19.11 (2009): 1163-1175.

[13] Gutenkunst, Ryan N., et al. "Universally sloppy parameter sensitivities in systems biology models." PLoS computational biology 3.10 (2007): e189.

[14] Machta, Benjamin B., et al. "Parameter space compression underlies emergent theories and predictive models." Science 342.6158 (2013): 604-607.

[15] Siebert, Stefan, et al. "Stem cell differentiation trajectories in Hydra resolved at single-cell resolution." Science 365.6451 (2019).

[16] C.T. Fincher et al, Cell type transcriptome atlas for the planarian Schmidtea mediterranea, Science. 2018 May 25;360(6391).

[17] Han, Xiaoping, et al. "Construction of a human cell landscape at single-cell level." Nature 581.7808 (2020): 303-309.

[18] Junqueira, Carneiro, Kelley, Basic Histology, 7th edition 1992 Appleton and Lange, p.66

https://doi.org/10.7554/eLife.76757.sa2

Article and author information

Author details

Alkesh Yadav

Raman Research Institute, Bangalore, India

Contribution
Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing - original draft, Writing - review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4268-8873
Quentin Vagne

Laboratoire Physico Chimie Curie, Institut Curie, CNRS UMR168, Paris, France

Contribution
Investigation, Methodology, Writing - review and editing

Competing interests
No competing interests declared
Pierre Sens

Laboratoire Physico Chimie Curie, Institut Curie, CNRS UMR168, Paris, France

Contribution
Funding acquisition, Project administration, Writing - review and editing

Competing interests
Reviewing editor, eLife

"This ORCID iD identifies the author of this article:" 0000-0003-4523-3791
Garud Iyengar

Industrial Engineering and Operations Research, Columbia University, New York, United States

Contribution
Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing - original draft, Writing - review and editing

For correspondence
garud@ieor.columbia.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-6546-4154
Madan Rao

Simons Centre for the Study of Living Machines, National Centre for Biological Sciences, Bangalore, India

Contribution
Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing - original draft, Writing - review and editing

For correspondence
madan@ncbs.res.in

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-6210-6386

Funding

Department of Atomic Energy, Government of India (RTI4006)

Madan Rao

Simons Foundation (287975)

Madan Rao

JC Bose Fellowship (DST-SERB)

Madan Rao

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank M Thattai, A Jaiman, S Ramaswamy, and A Varki for discussions, and S Krishna and R Bhat for very useful suggestions on the manuscript. We thank our group members at the Simons Centre for many incisive inputs. We are very grateful to P Babu and PS Sabarinath for consultations on the MSMS data and literature. MR acknowledges support from the Department of Atomic Energy (India), under project no. RTI4006, the Simons Foundation (grant no. 287975), and a JC Bose Fellowship from DST-SERB (India). We acknowledge the computational facilities at NCBS. This work has received support under the program Investissements d’Avenir launched by the French Government and implemented by ANR with the references ANR-10-LABX-0038 and ANR-10-IDEX-0001-02 PSL. MR thanks Institut Curie for hosting a visit under the Labex program, and QV thanks the Simons Centre (NCBS) for hosting his visit.

Senior Editor

Suzanne R Pfeffer, Stanford University School of Medicine, United States

Reviewing Editor

Felix Campelo, The Barcelona Institute of Science and Technology, Spain

Publication history

Preprint posted: May 19, 2020 (view preprint)
Received: January 4, 2022
Accepted: January 31, 2022
Accepted Manuscript published: February 17, 2022 (version 1)
Version of Record published: May 31, 2022 (version 2)

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

782

Page views
210

Downloads
2

Citations

Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Article PDF

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Alkesh Yadav
Quentin Vagne
Pierre Sens
Garud Iyengar
Madan Rao

(2022)

Glycan processing in the Golgi as optimal information coding that constrains cisternal number and enzyme specificity

eLife 11:e76757.

https://doi.org/10.7554/eLife.76757

Share this article

Cite this article

Living cells display a complex glycan distribution.

Enzymatic reaction and transport network in the secretory pathway.

Trade-offs amongst the glycan synthesis parameters, enzyme specificity σ, cisternal number NC, and enzyme number NE to achieve a complex target distribution c∗.

Fidelity of glycan distribution and optimal enzyme properties to achieve a complex target distribution.

Optimal enzyme partitioning in cisternae.

Stiff and sloppy directions in the optimization parameters.

Strategies for achieving high glycan diversity.

The binned mass spectroscopy (MS) data (blue) approximates the raw MS data (red) very well.

Log likelihood vs. number of components (N) in the Gaussian mixture model (GMM).

Enzyme parameters taken from Table 3 in Umaña and Bailey, 1997 that we use to calculate the bounds on the reaction rate R.

Flow chart showing the optimization schemes for Optimization A and B.

Glycan concentration profile calculated from the model using (a) formula (31) for NE=NC=1 and (b) formulae (32)–(36) for NE=1,NC=2.

Glycan profile {ck:k=1,…,Ns} as a function of specificity σ (a, c) and reaction rates R (b, d).

Optimum fidelity D¯K⁢L as a function of (NE,NC) for different values of ℓb/Ns, where ℓb bounds the deformation in the ideal length ℓα(0) of an enzyme α=1,…,NE.

Recovering the σ values for different target distribution.

Distribution of local minima.

D¯ for various initial conditions, sorted in increasing order for clarity.

Diversity vs. NC for different values of σ keeping NE=1 fixed, for three different values of the threshold, ct⁢h=1Ns,12⁢Ns,14⁢Ns.

Table of symbols and their definitions.

Author details

Alkesh Yadav

Contribution

Competing interests

Quentin Vagne

Contribution

Competing interests

Pierre Sens

Contribution

Competing interests

Garud Iyengar

Contribution

For correspondence

Competing interests

Madan Rao

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Further reading

Trade-offs amongst the glycan synthesis parameters, enzyme specificity $σ$ , cisternal number $N_{C}$ , and enzyme number $N_{E}$ to achieve a complex target distribution $c^{*}$ .

Log likelihood vs. number of components ( $N$ ) in the Gaussian mixture model (GMM).

Enzyme parameters taken from Table 3 in Umaña and Bailey, 1997 that we use to calculate the bounds on the reaction rate $R$ .

Glycan concentration profile calculated from the model using (a) formula (31) for $N_{E} = N_{C} = 1$ and (b) formulae (32)–(36) for $N_{E} = 1, N_{C} = 2$ .

Glycan profile ${c_{k} : k = 1, \dots, N_{s}}$ as a function of specificity $σ$ (a, c) and reaction rates $R$ (b, d).

Optimum fidelity ${\bar{D}}_{K L}$ as a function of $(N_{E}, N_{C})$ for different values of $ℓ_{b} / N_{s}$ , where $ℓ_{b}$ bounds the deformation in the ideal length $ℓ_{α}^{(0)}$ of an enzyme $α = 1, \dots, N_{E}$ .

Recovering the $σ$ values for different target distribution.

$\bar{D}$ for various initial conditions, sorted in increasing order for clarity.

Diversity vs. $N_{C}$ for different values of $σ$ keeping $N_{E} = 1$ fixed, for three different values of the threshold, $c_{t h} = \frac{1}{N_{s}}, \frac{1}{2 N_{s}}, \frac{1}{4 N_{s}}$ .