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

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)

Figure 2

Download asset Open asset

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.

We consider an array of $N_C$ Golgi cisternae, labelled by $j=1,\mathrm{\dots },N_C$, between the ER and PM (Figure 2). GBPs, denoted as ${\displaystyle \mathcal{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 ${\displaystyle \mathcal{P}{c}_1^{(1)}}$ reacts with ${\displaystyle c_0^{(1)}}$ to form the glycosylated acceptor ${\displaystyle \mathcal{P}{c}_2^{(1)}}$, following an MM reaction (1) catalyzed by the appropriate enzyme. The acceptor ${\displaystyle \mathcal{P}{c}_2^{(1)}}$ has the potential of being transformed into ${\displaystyle \mathcal{P}{c}_3^{(1)}}$, and so on, provided the requisite enzymes are present in that cisterna. This leads to the sequence of enzymatic reactions ${\displaystyle \mathcal{P}{c}_1^{(1)}\to \mathcal{P}{c}_2^{(1)}\to \dots \mathcal{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 ${\displaystyle {\mu }^{(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_{\alpha }^{(j)}$, with $\alpha =1,\mathrm{\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\gg 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$, ${\displaystyle N_E}$ = 10–20, and ${\displaystyle 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_{\alpha }^{(j)}$ on the substrate ${\displaystyle \mathcal{P}{c}_k^{(j)}}$ in cisterna $j$ is given by

where $k=1,\mathrm{\dots }{N}_s-1$. In general, the forward, backward, and catalytic rates ${\omega }_f$, ${\omega }_b$, and ${\omega }_c$, respectively, depend on the cisternal label $j$, the reaction label $k$, and the enzyme label $\alpha$, 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 ${\omega }_f$ depends on the binding energy of the enzyme $E_{\alpha }^{(j)}$ to acceptor substrate ${\displaystyle \mathcal{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 ${\omega }_f(j,k,\alpha )\propto P^{(j)}(k,\alpha )$, where $P^{(j)}(k,\alpha )\in (0,1)$ is the binding probability of enzyme $E_{\alpha }^{(j)}$ with substrate ${\displaystyle \mathcal{P}{c}_k^{(j)}}$, and define the binding probability $P^{(j)}(k,\alpha )$ 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 ${\mathbf{\ell }}_{\alpha }^{(j)}$ and ${\mathbf{\ell }}_k$ denote, respectively, the optimal ‘shape’ for enzyme $E_{\alpha }^{(j)}$ and the substrate ${\displaystyle \mathcal{P}{c}_k^{(j)}}$. We assume that the mismatch (or distortion) energy between the substrate $k$ and enzyme $E_{\alpha }^{(j)}$ is $\parallel {\mathbf{\ell }}_k-{\mathbf{\ell }}_{\alpha }^{(j)}\parallel$, with a binding probability given by

where $\parallel \cdot \parallel$ is a distance metric defined on the space of ${\mathbf{\ell }}_{\alpha }^{(j)}$ (e.g. the square of the ${\mathrm{\ell }}_2$-norm would be related to an elastic distortion model [Savir and Tlusty, 2007]) and the vector ${\displaystyle \mathit{\sigma }\equiv [{\sigma }_{\alpha }^{(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 ${\sigma }_{\alpha }^{(j)}$ indicates a highly specific enzyme, a small value of ${\sigma }_{\alpha }^{(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

where $c_k^{(j)}$ is the concentration of the acceptor substrate ${\displaystyle \mathcal{P}{c}_k^{(j)}}$ in compartment $j$.

Together with the constancy of the total enzyme concentration, ${\displaystyle {\left[E_{\alpha }^{(j)}\right]}_{tot}=E_{\alpha }^{(j)}+\sum _{k=1}^{N_s}}$ ${\displaystyle \left[E_{\alpha }^{(j)}-\mathcal{P}{c}_k^{(j)}-c_0^{(j)}\right]}$, this immediately fixes the kinetics of product formation (not including inter-cisternal transport),

where

and

This reparametrization of the reaction rates ${\omega }_f,{\omega }_b,{\omega }_c$ in terms of ${\displaystyle \mathbf{M},\mathbf{V}}$ is convenient since it relates to experimentally measurable parameters $V_{max}$ and MM constant $K_M$, for each $(j,k,\alpha )$, which can be easily read out (see Appendix 2). As is the usual case, the maximum velocity $V_{max}$ is not an intrinsic property of the enzyme because it is dependent on the enzyme concentration ${\left[E_{\alpha }^{(j)}\right]}_{tot}$; while $K_M(j,k,\alpha )=M(j,k,\alpha )c_0^{(j)}/P^{(j)}(k,\alpha )$ is an intrinsic parameter of the enzyme and the enzyme-substrate interaction. The enzyme catalytic efficiency, the so-called “${\displaystyle k_{cat}/K_M}$” ${\displaystyle \propto P^{(j)}(k,\alpha )}$, 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 ${\mu }^{(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, ${\displaystyle \mathbf{c}\equiv c_k^{(j)}}$, as a function of the independent vectors ${\displaystyle \mathbf{M}\equiv [M(j,k,\alpha )]}$, ${\displaystyle \mathbf{V}\equiv [V(j,k,\alpha )]}$, and ${\displaystyle \mathbf{L}\equiv [P^{(j)}(k,\alpha )]}$, the transport rates ${\displaystyle \mathit{\mu }\equiv [{\mu }^{(j)}]}$ and specificity, ${\displaystyle \mathit{\sigma }\equiv [{\sigma }_{\alpha }^{(j)}]}$.

Let ${\displaystyle {\mathit{c}}^{\ast }}$ 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 ${\displaystyle {\mathit{c}}^{\ast }}$. Let ${\displaystyle \overline{\mathit{c}}}$ denote the normalized steady-state glycan concentration distribution displayed on the PM. Then Equation 16 implies that ${\overline{c}}_k={\mu }^{(N_C)}{c}_k^{(N_C)}$, $k=1,\mathrm{\dots },N_s$. We measure the fidelity ${\displaystyle F({\mathit{c}}^{\ast }\Vert \overline{\mathit{c}})}$ between the ${\displaystyle {\mathbf{c}}^{\mathbf{\ast }}}$ and ${\displaystyle \overline{\mathbf{c}}}$ by the ratio of the KL divergence ${\displaystyle D({\mathit{c}}^{\ast }\Vert \mathit{c})}$ (Cover and Thomas, 2012; MacKay, 2003) to the entropy ${\displaystyle H({\mathit{c}}^{\ast })}$

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 ${\displaystyle F({\mathit{c}}^{\ast }\Vert \overline{\mathit{c}})}$, vice versa.

Thus, the problem of designing a sequential synthesis mechanism that approximates ${\displaystyle {\mathbf{c}}^{\ast }}$ for a given enzyme specificity ${\displaystyle \mathit{\sigma }}$, transport rate ${\displaystyle \mathit{\mu }}$, number of enzymes $N_E$, and number of cisternae $N_C$ is given by

where we emphasize that the optimum fidelity ${\displaystyle \overline{D}(\sigma ,N_E,N_C,{\mathit{c}}^{\mathbf{\ast }})}$ is a function of ${\displaystyle (\sigma ,N_E,N_C,{\mathit{c}}^{\mathbf{\ast }})}$. 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 ${\displaystyle (\mathit{\mu },\mathbf{M},\mathbf{V},\mathbf{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 ${\displaystyle \mathit{\mu },\mathbf{R}}$, and $\mathbf{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.

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 ${\mathbf{T}}^{(i)}$ denote the $i$-component GMM approximation for the human T-cell MSMS data. We sample this target distribution at indices $k=1,\mathrm{\dots },N_s$, that represent the glycan indices, and then renormalize to obtain the discrete distribution $\{T_k^{(i)},k=1,\mathrm{\dots },N_s\}$. To highlight the role of target distribution complexity, we focus on the 3-GMM ${\mathit{T}}^{(3)}$ (low complexity) and 20-GMM approximation ${\mathit{T}}^{(20)}$ (high complexity) in describing our results.

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

Figure 4

Download asset Open asset

1. The optimal fidelity $\overline{D}(\sigma ,N_E,N_C,{\mathbf{c}}^{\ast })$ is a convex function of $\sigma$ for fixed values for other parameters (see Figure 3), that is, it first decreases with $\sigma$ and then increases beyond a critical value of ${\sigma }_{\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 $\overline{D}(\sigma ,N_E,N_C,{\mathbf{c}}^{\ast })$ is decreasing in $N_C$ and $N_E$ for fixed values of the other parameters, and increasing in the complexity of ${\mathbf{c}}^{\ast }$ for fixed $(\sigma ,N_C)$. The marginal contribution of $N_C$ and $N_E$ in improving fidelity $\overline{D}$ is approximately equal (see Figure 4a and b). We discuss the origin of this symmetry later in this section.

2. The optimal enzyme specificity ${\sigma }_{\min }({\mathbf{c}}^{\ast },N_C)={\text{argmin}}_{\sigma }\{\overline{D}(\sigma ,{\overline{N}}_E,N_C,{\mathbf{c}}^{\ast })\}$, which minimizes the error as function of $(N_C,{\mathbf{c}}^{\ast })$ with $N_E$ fixed at ${\overline{N}}_E$, is an increasing function of $N_C$ and the complexity of the target distribution ${\mathbf{c}}^{\ast }$ (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 $\sigma$: since a GMM approximation with fewer peaks has wider peaks, ${\sigma }_{\min }$ is low, and vice versa. Similar results hold when $N_C$ is fixed at ${\overline{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 $\overline{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.

Having studied the optimum $N_E,N_C,\sigma$ 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 $\sigma$ 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 $\mathbf{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_{eff}(j,k)$ for $P{c}_k\to P{c}_{k+1}$ in the jth cisterna as

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 ${\displaystyle {\mathbf{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

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).

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 $\mathbf{R}$, $\mathit{\mu }$, $\mathbf{L}$, and $\sigma$ is typically rugged. We analyse the geometry of this fitness landscape by doing a local Hessian analysis about the optimal solutions.

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 ($\mathbf{R}$)-transport rate ($\mathit{\mu }$) 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.

We analyse the change in fidelity on small perturbations in $\mathbf{R}$, $\mathit{\mu }$, $\mathbf{L}$, and $\sigma$ 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 $\sigma$ change the glycan profile a lot more compared to perturbations in the other parameters and this stiffness in $\sigma$ generally decreases on increasing $N_E,N_C$ (Figure 6a–c). Small perturbations in $\mathit{\mu }$ and some $\mathit{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 $\sigma$. The glycan profile is robust to perturbations in most $\mathit{R}$ and some $\mathit{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

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.

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_{th}$. This last condition is necessary because very small peaks will not be distinguishable in the presence of noise. In computing the diversity from our chemical synthesis model, we have chosen the threshold to be $c_{th}=1/N_s$, where $N_s$ is the total number of glycan species. We have checked that the qualitative results do not depend on this choice (see Appendix 10—figure 1).

We use the sigmoid function ${(1+e^{-x/\tau })}^{-1}$ as a differentiable approximation to the Heaviside function $\mathrm{\Theta }(x)$ and define the following optimization to maximize diversity for a given set of parameter values, $N_E,N_C,\sigma$:

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

The results displayed in Figure 7a show that for a fixed specificity $\sigma$ the diversity at first increases with the number of cisternae $N_C$, and then saturates at a value that depends on $\sigma$. 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

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.

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 ${\overline{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:

1. 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 }{\overline{I}}_k$. The plots in Appendix 1—figure 2a show that the results are not sensitive to the specific value of $I_{\max }$.

2. 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_{max}}}=1-e^{-{\overline{I}}_k}$

   which adequately suppresses this chemical noise.

3. Bootstrapped glycan data: The count n_k at the glycan index $k$ is distributed according to the following distribution:

   $n_k=\{\begin{array}{ccc}0\hfill & (1-p_k)\hfill & n=0\hfill \\ n\hfill & p_k{e}^{-I_k}{\displaystyle \frac{{(I_k)}^n}{n!}}\hfill & n\ge 1.\hfill \end{array}$

   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

Appendix 1—figure 2

Download asset Open asset

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_{GMM}(\theta )={\sum }_{i=1}^m{w}_i{N}_{{\eta }_i,{\mathrm{\Delta }}_i}$, where $N_{{\eta }_i,{\mathrm{\Delta }}_i}$ denotes the Gaussian distribution with mean ${\eta }_i$ and variance ${\mathrm{\Delta }}_i$, $w_i\ge 0$ and ${\sum }_{i=1}^m{w}_i=1$, the parameter vector $\mathit{\theta }=(\mathit{w},\mathit{\eta },\mathbf{\Delta })$ . We compute the optimal $m$-component GMM approximation by minimizing the KL divergence $D(p_T||p_{GMM}(\mathit{\theta }))$ as a function of parameter vector $\mathit{\theta }$. Since

the optimization problem ${\min }_{\mathit{\theta }}D(p_T||p_{GMM}(\mathit{\theta }))$ is equivalent to

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

Then $z_i^{(t)}(k)\ge 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

Then, we have that

and

Define

Then, we have that

Therefore, the iterative algorithm in (Equation 11) generates a sequence $\{{\mathit{\theta }}^{(t)}:t\ge 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.

1. $w$-update

   (9) ${\mathit{w}}^{(t+1)}={\displaystyle \text{arg}\underset{\mathit{w}}{max}\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)}}$

2. $\eta$-update

   (10) ${\eta }_i^{(t+1)}={\displaystyle \text{arg}\underset{\eta }{min}\sum _{k=1}^{N_S}{p}_T(k){\hat{z}}_i(k)|k-{\eta }_i{|}^2⟹{\eta }_i^{(t+1)}=\frac{\sum _{k=1}^{N_s}{z}_i^{(t)}(k)k}{\sum _{k=1}^{N_s}{z}_i^{(t)}(k)}.}$

3. $\mathrm{\Delta }$-update

   (11) $\begin{array}{ll}{\mathrm{\Delta }}_i^{(t+1)}& ={\displaystyle {\text{argmax}}_{\mathrm{\Delta }\ge {\mathrm{\Delta }}_{cut}}\left\{-\frac{\sum _{k=1}^{N_s}{p}_T(k)z_i^{(t)}(k)|k-{\eta }_i^{t+1}{|}^2}{2\mathrm{\Delta }}-\log (\mathrm{\Delta })\right\}}\\ \\ & ={\displaystyle max\left(\sqrt{\frac{\sum _{k=1}^{N_S}{p}_T(k)z_i^{(t)}(k)|k-{\eta }_i^{(t+1)}{|}^2}{\sum _{k=1}^{N_s}{p}_T(k)z_i^{(t)}(k)}},{\mathrm{\Delta }}_{cut}\right),}\end{array}$

where ${\mathrm{\Delta }}_{cut}$ is the minimum allowed width of the Gaussians, in our case ${\mathrm{\Delta }}_{cut}=1$ since glycan index, $k\in \{1,2,\mathrm{\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_{GMM})$), 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.

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 $\mu =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

where $K_M^{(\alpha )}$ and $V_{\max }^{(\alpha )}$ 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 ($\nu$) to be 2.5 µm³ (Umaña and Bailey, 1997; Krambeck et al., 2009). This gives

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

On including the rates of injection ($q$) and inter-cisternal transport ${\mu }^{(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

for cisterna-1, and

for cisternae $j=2,3,\mathrm{\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\ge 0$. Equations (13)-(14) automatically obey the conservation law for the protein concentration ($p$), that is, the total protein concentration $p^{(j)}={\sum }_{k^{\prime }=1}^{N_s}{c}_{k^{\prime }}^{(j)}$ in the jth cisterna automatically satisfies

for $j=2,3,\mathrm{\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,\alpha )=V(j,k,\alpha )/q$ and ${\mu }^{(j)}={\mu }^{(j)}/q$, gives the following recursion relations for the steady-state concentrations of the glycans in each cisterna. In the first cisterna,