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

Alkesh Yadav  Raman Research Institute, Bangalore 560080, India  Quentin Vagne
Laboratoire Physico Chimie Curie, Institut Curie, CNRS UMR168, 75005 Paris,
France  Pierre Sens  Laboratoire Physico Chimie Curie, Institut Curie, CNRS
UMR168, 75005 Paris, France  Garud Iyengar  Industrial Engineering and
Operations Research, Columbia University, New York 10027, USA  Madan Rao
Simons Centre for the Study of Living Machines, National Centre for Biological
Sciences (TIFR), Bangalore 560065, India<EMAIL_ADDRESS>
###### 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
realises the desired target glycan distribution for a particular cell type and
niche. In this paper, we quantitatively analyse the tradeoffs between the
number of cisternae and the number and specificity of enzymes, in order to
synthesize a prescribed target glycan distribution of a certain complexity. We
find that to synthesize complex distributions, such as those observed in real
cells, one needs to have multiple cisternae and precise enzyme partitioning in
the Golgi. Additionally, for fixed number of enzymes and cisternae, there is
an optimal level of specificity of enzymes that achieves the target
distribution with high fidelity. Our results show how the complexity of the
target glycan distribution places functional constraints on the Golgi
cisternal number and enzyme specificity.

##  I Introduction

A majority of the proteins synthesized in the endoplasmic reticulum (ER) are
transferred to the Golgi cisternae for further chemical modification by
glycosylation Alberts  _et al._ (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  _et al._ (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  _et al._
(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, 1998);
Pothukuchi  _et al._ (2019).

In this paper, we will focus on the role of glycans as markers of cell
identity. For the glycans to play this role, they must inevitably represent a
molecular code 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). In this paper, we study one
aspect of molecular coding, namely the fidelity of this molecular code
generation. While it has been recognised that fidelity of the glycan code is
necessary for reliable cellular recognition Demetriou  _et al._ (2001), a
quantitative measure of fidelity of the code and its implications for cellular
structure and organization are lacking.

There are two aspects of the cell-type specific glycan code that have an
important bearing on quantifying fidelity. The first is that extant glycan
distributions have high complexity, owing to evolutionary pressures arising
from (a) reliable cell type identification amongst a large set of different
cell types in a complex organism, the preservation and diversification of
“self-recognition” Drickamer and Taylor (1998), (b) pathogen-mediated
selection pressures Varki  _et al._ (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). Here, we will interpret this to mean that the target
distribution of glycans of a given cell type is complex; in Sect. II we define
a quantitative measure for complexity and demonstrate its implications in the
context of human T-cells. The second is that the cellular machinery for the
synthesis of glycans, which involves sequential chemical processing via
cisternal resident enzymes and cisternal transport, is subject to variation
and noise Varki (2017, 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 Sect. III.

In this paper, 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. This KL divergence
is a function of the cellular parameters governing glycan synthesis, such as
the number and specificity of enzymes, inter-cisternal transfer rates, and
number of cisternae (Sect. V). We analyze the tradeoffs between the number of
cisternae and the number and specificity of enzymes, in order to achieve a
prescribed target glycan distribution with high fidelity (Sect. VI). Our
analysis leads to a number of interesting results, of which we list a few
here: (i) In order to construct an accurate representation of a complex target
distribution, such as those observed in real cells, one needs to have multiple
cisternae and precise enzyme partitioning. Low complexity target distributions
can be achieved with fewer cisternae. (ii) This definition of fidelity of the
glycan code, allows us to provide a quantitative argument for the evolutionary
requirement of multiple-compartments. (iii) For fixed number of enzymes and
cisternae, there is an optimal level of specificity of enzymes that achieves
the complex target distribution with high fidelity. (iv) Keeping the number of
enzymes fixed, having low specificity or sloppy enzymes and larger cisternal
number could give rise to a diverse repertoire of functional glycans, a
strategy used in organisms such as plants and algae.

Stated another way, our results imply that the pressure to achieve the target
glycan code for a given cell type, places strong constraints on the cisternal
number and enzyme specificity Sengupta and Linstedt (2011). This would suggest
that a description of the nonequilibrium assembly of a fixed number of Golgi
cisternae must combine the dynamics of chemical processing with membrane
dynamics involving fission, fusion and transport Sachdeva  _et al._ (2016);
Sens and Rao (2013), opening up a new direction for future research.

##  II Complexity of glycan code in real cells

(a)

(b)

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

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. A more complex or richer class
of glycan profiles is able to support a larger number of well separated
profiles, and therefore, a larger number cell types, or equivalently, a more
complex organism111A rigorous definition of complexity can be given in terms
of the Kullback-Leibler metric Cover and Thomas (2012); MacKay (2003) between
two glycan profiles. We declare that two profiles are distinguishable only if
the Kullback-Leibler 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 Kullback-Leibler distance of any pair of profiles
is larger than the tolerance.

In order to implement a quantitative measure of complexity, we first need a
consistent way of smoothening or coarse-graining the discrete glycan
distribution to remove measurement and synthesis noise. In this paper, we
approximate the glycan profile as mixture of Gaussian densities with specified
number of components that are supported on a finite set of indices Bacharoglou
(2010). Since the complexity of $k$-component Gaussian is an increasing
function of $k$, we use the number of component $k$ and complexity
interchangeably.

Using this definition we demonstrate that the glycan profiles of typical
mammalian cells are very complex. We obtain target profiles for a given cell
type from Mass Spectrometry coupled with determination of molecular structure
(MSMS) measurements Cummings and Crocker (2020). Fig. 1 shows the the MSMS
data from human T-cells and human and mouse neutrophils Cummings and Crocker
(2020), and their coarse-grained representations using Gaussian mixture models
(GMM) of differing complexity - a low complexity $k=3$ GMM and high complexity
$k=20$ GMM. It is clear from Fig. 1, that the more complex $k=20$ GMM is a
better representation of the MSMS data as compared to the less complex $k=3$
GMM. Indeed the $k=20$ Gaussian mixture model is the best compromise between
faithfulness of the representation and cost of an additional component, as
seen from the saturation of the likelihood function MacKay (2003). Details of
this systematic coarse-graining procedure appear in Sect. VI.2 and Appendix G.

Having demonstrated the complexity of the typical glycan distributions
associated with a given cell type, we will now describe a general model of the
cellular machinery that is capable of synthesizing glycans of the desired
complexity. We expect that cells need a more elaborate mechanism to produce
profiles from a more complex set.

##  III Synthesis of glycans in the Golgi cisternae

The glycan display at the cell surface is a result of proteins that flux
through and undergo sequential chemical modification in the secretory pathway,
comprising an array of Golgi cisternae situated between the ER and the PM, as
depicted in Fig. 2. Glycan-binding proteins (GBPs) are delivered from the ER
to the first cisterna, whereupon they are processed by the resident enzymes in
a sequence of steps that constitute the N-glycosylation process Varki  _et
al._ (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, e.g. UDP-Gal, is transferred to the acceptor
substrate, by a Michaelis-Menten (MM) type reaction Varki  _et al._ (2009)

$\displaystyle\small{\mbox{Acceptor}+\mbox{glycosyl donor}+\mbox{Enzyme}}$ (1)
$\displaystyle\small{\xrightleftharpoons[\omega_{b}]{\omega_{f}}}\quad\small{\left[\mbox{Acceptor}\cdot\mbox{glycosyl
donor}\cdot\mbox{Enzyme}\right]}\quad\small{\xrightarrow[]{\omega_{c}}}\quad\small{\mbox{glycosylated
acceptor}+\mbox{nucleotide}+\mbox{Enzyme}}$

From the first cisterna, the proteins with attached sugars are delivered to
the second cisterna at a given inter-cisternal transfer rate, where further
chemical processing catalysed by the enzymes resident in the second cisterna
occurs. This chemical processing and inter-cisternal transfer continues until
the last cisterna, thereupon the fully processed glycans are displayed at the
PM Varki  _et al._ (2009). The network of chemical processing and inter-
cisternal transfer forms the basis the physical model that we will describe
next.

Any physical model of such a network of enzymatic reactions and cisternal
transfer needs to be augmented by reaction and transfer rates and chemical
abundances. To obtain the range of allowable values for the reaction rates and
chemical abundances, we use the elaborate enzymatic reaction models, such as
the KB2005 model 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, and compare these predictions with
N-glycan mass spectrum data. For the allowable rates of cisternal transfer, we
rely on the recent study by Ungar and coworkers 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.

Figure 2: Enzymatic reaction and transport network in the secretory pathway.
Represented here is the array of Golgi cisternae (blue) indexed by
$j=1,\ldots,N_{C}$ situated between the ER and PM. Glycan-binding proteins
${\cal{P}}c_{1}^{(1)}$ are injected from the ER to cisterna-1 at rate $q$.
Superimposed is transition network of chemical reactions (column) -
intercisternal transfer (rows), the latter with rates $\mu^{(j)}$.
${\cal{P}}c^{(j)}_{k}$ denotes the acceptor substrate in compartment $j$ and
the glycosyl donor $c_{0}$ is chemostated in each cisterna. This results in a
frequency distribution of glycans displayed at the PM (red curve), that is
representative of the cell type.

##  IV Model Definition

###  IV.1 Chemical reaction and transport network in cisternae

With this background, we now define our quantitative model for chemical
processing and transport in the secretory pathway. We consider an array of
$N_{C}$ Golgi cisternae, labelled by $j=1,\ldots,N_{C}$, between the ER and
the PM (Fig. 2). Glycan-binding proteins (GBPs), denoted as
${\cal{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 $j$-th cisterna is chemostated Varki  _et al._ (2009);
Hirschberg  _et al._ (1998); Caffaro and Hirschberg (2006); Berninsone and
Hirschberg (2000), thus in our model we hold its concentration $c^{(j)}_{0}$
constant in time for each $j$. The acceptor ${\cal{P}}c_{1}^{(1)}$ reacts with
$c_{0}^{(1)}$ to form the glycosylated acceptor ${\cal{P}}c_{2}^{(1)}$,
following an MM-reaction (1) catalysed by the appropriate enzyme. The acceptor
${\cal{P}}c_{2}^{(1)}$ has the potential of being transformed into
${\cal{P}}c_{3}^{(1)}$, and so on, provided the requisite enzymes are present
in that cisterna. This leads to the sequence of enzymatic reactions
${\cal{P}}c_{1}^{(1)}\to{\cal{P}}c_{2}^{(1)}\to\ldots{\cal{P}}c_{k}^{(1)}\to\ldots$,
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 $\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 Fig. 2.
Although, in this paper, we focus on a sequence of reactions that form a line-
graph, the methodology we propose extends to tree-like reaction sequences, and
more generally to reaction sequences that form a directed acyclic graphs
Trinajstic (2018).

Let $N_{s}$ denote the maximum number of possible glycosylation reactions in
each cisterna $j$, catalysed by enzymes labelled as $E^{(j)}_{\alpha}$, with
$\alpha=1,\ldots,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 Fig. 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 localisation, by setting their
concentrations to zero in those cisternae where they are not present.

Now the action of enzyme $E^{(j)}_{\alpha}$ on the substrate
${\cal{P}}c^{(j)}_{k}$ in cisterna $j$ is given by

$\begin{split}{\cal{P}}c_{k}^{(j)}+E_{\alpha}^{(j)}\xrightleftharpoons[\omega_{b}(j,k,\alpha)]{\omega_{f}(j,k,\alpha)c_{0}^{(j)}}\left[E_{\alpha}^{(j)}-{\cal{P}}c_{k}^{(j)}-c_{0}^{(j)}\right]&\xrightarrow{\omega_{c}(j,k,\alpha)}{\cal{P}}c_{k+1}^{(j)}+E_{\alpha}^{(j)}\end{split}$
(2)

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$, that parametrise
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}$ to depend on the
binding energy of the enzyme $E_{\alpha}^{(j)}$ to acceptor substrate
${\cal{P}}c_{k}^{(j)}$ and a physical variable that characterises
cisterna-$j$.

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 their
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^{(j)}_{\alpha}$ with substrate ${\cal{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), of enzyme-substrate induced fit.

Let $\bm{l}^{(j)}_{\alpha}$ and $\bm{l}_{k}$ denote, respectively, the optimal
“shape” for enzyme $E^{(j)}_{\alpha}$ and the substrate
${\cal{P}}c_{k}^{(j)}$. We assume that the mismatch (or distortion) energy
between the substrate $k$ and enzyme $E^{(j)}_{\alpha}$ is
$\|\bm{l}_{k}-\bm{l}_{\alpha}^{(j)}\|$, with a binding probability given by,

$P^{(j)}(k,\alpha)=\exp\left(-\sigma_{\alpha}^{(j)}\|\bm{l}_{k}-\bm{l}_{\alpha}^{(j)}\|\right)$
(3)

where $\|.\|$ is a distance metric defined on the space of
$\bm{l}_{\alpha}^{(j)}$ (e.g., the square of the $\ell_{2}$-norm would be
related to an elastic distortion model Savir and Tlusty (2007)) and the vector
$\bm{\sigma}\equiv[\sigma_{\alpha}^{(j)}]$ parametrises enzyme specificity. A
large value of $\sigma_{\alpha}^{(j)}$ indicates a highly specific enzyme, a
small value of $\sigma_{\alpha}^{(j)}$ indicates a promiscuous or sloppy
enzyme. It is recognised that the degree of enzyme specificity or sloppiness
is an important determinant of glycan distribution Varki  _et al._ (2009);
Roseman (2001); Hossler  _et al._ (2007); Yang  _et al._ (2018).

As in Umaña and Bailey (1997); Krambeck  _et al._ (2009); Krambeck and
Betenbaugh (2005), our synthesis model is mean-field, in that we ignore
stochasticity in glycan synthesis that may arise from low copy numbers of
substrates and enzymes, multiple substrates competing for the same enzymes,
and kinetics of inter-cisternal transfer. Then the usual MM-steady state
condition on (2), which assumes that the concentration of the intermediate
enzyme-substrate complex does not change with time, implies

$\left[E_{\alpha}^{(j)}-{\cal{P}}c_{k}^{(j)}-c_{0}^{(j)}\right]=\frac{\omega_{f}(j,k,\alpha)\,c_{0}^{(j)}}{\omega_{b}(j,k,\alpha)+\omega_{c}(j,k,\alpha)}E_{\alpha}^{(j)}c_{k}^{(j)}.\textbf{{}}$

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

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

$\frac{dc_{k+1}^{(j)}}{dt}=\sum_{\alpha=1}^{N_{E}}\frac{V(j,k,\alpha)P^{(j)}(k,\alpha)c_{k}^{(j)}}{M(j,k,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(j)}(k^{\prime},\alpha)c_{k^{\prime}}^{(j)}}{M(j,k^{\prime},\alpha)}\right)}$
(4)

where

$M(j,k,\alpha)=\frac{\omega_{b}(j,k,\alpha)+\omega_{c}(j,k,\alpha)}{\omega_{f}(j,k,\alpha)c_{0}^{(j)}}P^{(j)}(k,\alpha)$

and

$V(j,k,\alpha)=\omega_{c}(j,k,\alpha)\left[E_{\alpha}^{(j)}\right]_{tot}\,.$

From the above, the experimentally measurable parameters $V_{max}$ and MM-
constant $K_{M}$, for each $(j,k,\alpha)$ can be easily read out. As is the
usual case, the maximum velocity $V_{max}$ is not an intrinsic property of the
enzyme, because it is dependent on the enzyme concentration
$\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 $\text{``$k_{cat}/K_{M}$"}\propto
P^{(j)}(k,\alpha)$ and 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 A 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$). Rescaling the kinetic parameters in terms of the
injection rate $q$, i.e. $V(j,k,\alpha)={V(j,k,\alpha)}/{q}$ and
$\mu^{(j)}={\mu^{(j)}}/{q}$, we see that the steady state concentrations of
the glycans in each cisterna satisfy the following recursion relations (see,
Appendix A). In the first cisterna,

$\displaystyle c_{1}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{1}{\mu^{(1)}+\sum_{\alpha=1}^{N_{E}}\frac{V(1,1,\alpha)P^{(1)}(1,\alpha)c_{1}^{(1)}}{M(1,1,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(1)}(k^{\prime},\alpha)c_{k^{\prime}}^{(1)}}{M(1,k^{\prime},\alpha)}\right)}}$
$\displaystyle c_{k}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{\sum_{\alpha=1}^{N_{E}}\frac{V(1,k-1,\alpha)P^{(1)}(k-1,\alpha)c_{k-1}^{(1)}}{M(1,k-1,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(1)}(k^{\prime},\alpha)c_{k^{\prime}}^{(1)}}{M(1,k^{\prime},\alpha)}\right)}}{\mu^{(1)}+\sum_{\alpha=1}^{N_{E}}\frac{V(1,k,\alpha)P^{(1)}(k,\alpha)c_{k}^{(1)}}{M(1,k,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(1)}(k^{\prime},\alpha)c_{k^{\prime}}^{(1)}}{M(1,k^{\prime},\alpha)}\right)}}$
(5) $\displaystyle c_{N_{s}}^{(1)}$ $\displaystyle=$
$\displaystyle\frac{\sum_{\alpha=1}^{N_{E}}\frac{V(1,N_{s}-1,\alpha)P^{(1)}(N_{s}-1,\alpha)c_{N_{s}-1}^{(1)}}{M(1,N_{s}-1,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(1)}(k^{\prime},\alpha)c_{k^{\prime}}^{(1)}}{M(1,k^{\prime},\alpha)}\right)}}{\mu^{(1)}}$

and in cisternae $j\geq 2$,

$\displaystyle c_{1}^{(j)}$ $\displaystyle=$
$\displaystyle\frac{\mu^{(j-1)}c_{1}^{(j-1)}}{\mu^{(j)}+\sum_{\alpha=1}^{N_{E}}\frac{V(j,1,\alpha)P^{(j)}(1,\alpha)c_{1}^{(j)}}{M(j,1,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(j)}(k^{\prime},\alpha)c_{k^{\prime}}^{(j)}}{M(j,k^{\prime},\alpha)}\right)}}$
(6) $\displaystyle c_{k}^{(j)}$ $\displaystyle=$
$\displaystyle\frac{\mu^{(j-1)}c_{k}^{(j-1)}+\sum_{\alpha=1}^{N_{E}}\frac{V(j,k-1,\alpha)P^{(j)}(k-1,\alpha)c_{k-1}^{(j)}}{M(j,k-1,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(j)}(k^{\prime},\alpha)c_{k^{\prime}}^{(j)}}{M(j,k^{\prime},\alpha)}\right)}}{\mu^{(j)}+\sum_{\alpha=1}^{N_{E}}\frac{V(j,k,\alpha)P^{(j)}(k,\alpha)c_{k}^{(j)}}{M(j,k,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(j)}(k^{\prime},\alpha)c_{k^{\prime}}^{(j)}}{M(j,k^{\prime},\alpha)}\right)}}$
$\displaystyle c_{N_{s}}^{(j)}$ $\displaystyle=$
$\displaystyle\frac{\mu^{(j-1)}c_{N_{s}}^{(j-1)}+\sum_{\alpha=1}^{N_{E}}\frac{V(j,N_{s}-1,\alpha)P^{(j)}(N_{s}-1,\alpha)c_{N_{s}-1}^{(j)}}{M(j,N_{s}-1,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(j)}(k^{\prime},\alpha)c_{k^{\prime}}^{(j)}}{M(j,k^{\prime},\alpha)}\right)}}{\mu^{(j)}}$

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

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

These nonlinear recursion equations (5)-(6) have to be solved numerically to
obtain the steady state glycan concentrations, $\mathbf{c}\equiv c_{k}^{(j)}$,
as a function of the independent vectors $\mathbf{M}\equiv[M(j,k,\alpha)]$,
$\mathbf{V}\equiv[V(j,k,\alpha)]$, and $\mathbf{L}\equiv[P^{(j)}(k,\alpha)]$,
the transport rates $\bm{\mu}\equiv[\mu^{(j)}]$ and specificity,
$\bm{\sigma}\equiv[\sigma_{\alpha}^{(j)}]$.

##  V Optimization Problem

Now, with both the protocol for determining the target glycan distribution and
the sequential chemical processing model in hand, we can precisely define the
optimization problem referred to in the Introduction. Let $\mathbf{c}^{\ast}$
denote the “target” concentration distribution222We normalize the distribution
so that $\sum_{k=1}^{N_{s}}c^{\ast}_{k}=1$. for a particular cell type, i.e.
the goal of the sequential synthesis mechanism described in Sect. IV.1 is to
approximate $\mathbf{c}^{\ast}$. Let $\bar{\mathbf{c}}$ denote the steady
state glycan concentration distribution displayed on the PM - (6) implies that
$\bar{c}_{k}=\mu^{(N_{C})}c^{(N_{C})}_{k}$, $k=1,\ldots,N_{s}$. We measure the
fidelity between the $\bf{c}^{\ast}$ and $\bar{\bf{c}}$ by the Kullback-
Leibler metric Cover and Thomas (2012); MacKay (2003),

$D_{KL}(\mathbf{c}^{\ast}\|\bar{\mathbf{c}})=\sum_{k=1}^{N_{s}}c^{\ast}_{k}\ln\Big{(}\frac{c^{\ast}_{k}}{\bar{c}_{k}}\Big{)}=\sum_{k=1}^{N_{s}}c_{k}^{*}\ln\Big{(}\frac{c_{k}^{*}}{c_{k}^{(N_{C})}\mu^{(N_{C})}}\Big{)}$
(7)

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

$\text{\emph{Optimization A}}:\quad\min_{\mathbf{M},\ \mathbf{V},\ \mathbf{L}\
\geq\mathbf{0}}D_{KL}(\mathbf{c}^{\ast}\|\bar{\mathbf{c}})$ (8)

There is separation of time scales implicit in Optimization A – the chemical
kinetics of the production of glycans and their display on the PM happens over
cellular time scales, while the issues of tradeoffs and changes of parameters
are driven over evolutionary timescales.

Optimization A, though well-defined, is a hard problem, since the steady state
concentrations (6) are not explicitly known in terms of the parameters
$(\mathbf{M},\mathbf{V},\mathbf{L})$. In Appendix B, we formulate an
alternative problem _Optimization B_ in which the steady state concentrations
are defined explicitly in terms of a new parameters $\mathbf{R}$ and
$\mathbf{L}$, and in Appendix C we prove that Optimization A and Optimization
B are exactly equivalent. This is a crucial insight that allows us to obtain
all the results that follow.

In Appendix H, we describe the variant of the Sequential Quadratic Programming
(SQP) Boyd and Vandenberghe (2004), that we use to numerically solve the
optimization problem.

##  VI Results of optimization

To start with, the dimension of the optimization search space is extremely
large $\approx O(N_{s}\times N_{E}\times N_{C})$. To make the optimization
search more manageable, we ignore the $k$-dependence of the vectors
$(\mathbf{M},\mathbf{V})$, (or, alternatively of $\mathbf{R}$, see Appendix B
for details). The dependence on the reaction rates on the glycosyl substrate
is still present in the forward reactions via the enzyme-substrate binding
probability $P^{(j)}(k,\alpha)$. We further assume that shape function is a
number, $\bm{l}_{\alpha}^{(j)}=l_{\alpha}^{(j)}$ and that $\bm{l}_{k}=k$.
Finally we will drop the dependence of the specificity on $\alpha$ and $j$,
and take it to be a scalar $\sigma$. To fix our model, we will take the
distortion energy that appears in (3) to be the linear form
$|l_{k}-l_{\alpha}^{(j)}|$. Other metrics, such as
$|l_{k}-l_{\alpha}^{(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.

These restrictions significantly reduce the dimension of the optimization
search, so much so that in certain limits, we can solve the problem
analytically333In Appendix D we show that (19) can be solved analytically in
the limit $N_{s}\gg 1$, since the glycan index $k$ can be approximated by a
continuous variable, and the recursion relations for the steady state glycan
concentrations (5)-(6) can be cast as a matrix differential equation. This
allows us to obtain an _explicit_ expression for the steady state
concentration in terms of the parameters $(\mathbf{R},\mathbf{L})$. . This
helps us obtain some useful heuristics (Appendix E) on how to tune the
parameters, e.g. $N_{E}$, $N_{C}$, $\sigma$, and others, in order to generate
glycan distributions $\mathbf{c}$ of a given complexity. These heuristics
inform our more detailed optimization using “realistic” target distributions.

The calculations in Appendix E imply, as one might expect, that the synthesis
model needs to be more elaborate, i.e., needs a larger number of cisternae
$N_{C}$ or a larger number of enzymes $N_{E}$, in order to produce a more
complex glycan distribution. For a real cell type in a niche, the specific
elaboration of the synthesis machinery, would depend on a variety of control
costs associated with increasing $N_{E}$ and $N_{C}$. While an increase in the
number of enzymes would involve genetic and transcriptional costs, the costs
involved in increasing the number of cisternae could be rather subtle.

Notwithstanding the relative control costs of increasing $N_{E}$ and $N_{C}$,
it is clear from the special case, that increasing the number of cisternae
achieves the goal of obtaining an accurate representation of the target
distribution. Let us assume that the target distribution is
$c^{*}_{k}=\delta(k-M)$ for a fixed $M\gg 1$, i.e. $c^{\ast}_{k}=1$ when
$k=M$, and $0$ otherwise, and that the $N_{E}$ enzymes that catalyse the
reactions are highly specific. In this limit, Optimization A reduces to a
simple enumeration exercise Jaiman and Thattai (2018): clearly one needs
$N_{E}=M$, with one enzyme species for each of the $k=1,\ldots,M$ reactions,
in order to generate ${\cal{P}}c_{M}$. For a single Golgi cisterna with a
finite cisternal residence time (finite $\mu$), the chemical synthesis network
will generate a significant steady state concentration of lower index glycans
${\cal{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},\ldots,E_{M})$ with $E_{j}$ enzymes in cisterna
$j=1,\ldots,N_{c}$. This argument can be generalised to the case where the
target distribution is a finite sum of delta-functions. The more general case,
where the enzymes are allowed to have variable specificity, needs a more
detailed study, to which we turn to below.

###  VI.1 Target distribution from coarse-grained MSMS

As discussed in Sect. II, we obtain the target glycan distribution from glycan
profiles for real cells obtained using Mass Spectrometry coupled with
determination of molecular structure (MSMS) measurements 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 Sect. II, the Gaussian mixture models (GMM) 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 Gaussian mixture models (see Appendix G for details) to
create a hierarchy of increasingly complex distributions to approximate the
MSMS raw data. Thus, the $3$-GMM and $20$-GMM approximations represent the low
and high complexity benchmarks, respectively. In Appendix G, we show that the
likelihood for the glycan distribution of the human T-cell saturates at $20$
peaks. Thus, statistically speaking, the human T-cell glycan distribution is
accurately approximated by $20$ peaks.

This hierarchy allows us to study the trade-off between the complexity of the
target distribution and the complexity of the synthesis model needed to
generate the distribution as follows. Let $\mathbf{T}^{(i)}$ denote the
$i^{th}$-GMM approximation for the human T-cell MSMS data. We sample this
target distribution at indices $k=1,\ldots,N_{s}$, that represent the glycan
indices, and then renormalize to obtain the discrete distribution
$\\{T^{(i)}_{k},k=1,\ldots,N_{s}\\}$. Let
$H(\mathbf{T}^{(i)}):=-\sum_{k=1}^{N_{s}}T^{(i)}_{k}\log T^{(i)}_{k}$ denote
the entropy Cover and Thomas (2012) of the $i^{th}$-GMM approximation.
$H(\mathbf{T}^{(i)})$ quantifies statistical information in the target
distribution $\mathbf{T}^{(i)}$. We evaluate the fidelity of the distribution
generated by the synthesis model to this target distribution by the ratio of
the Kullback-Leibler distance to the entropy of the target distribution:

${\bar{D}}(\sigma,N_{E},N_{C},\mathbf{T}^{(i)}):=\frac{D(\sigma,N_{E},N_{C},\mathbf{T}^{(i)})}{H(\mathbf{T}^{(i)})}$
(9)

This normalization allows us evaluate the fidelity of the synthesis model to
the target distribution $\mathbf{T}^{(i)}$ as a fraction of the total
statistical information in the target distribution $\mathbf{T}^{(i)}$.

###  VI.2 Tradeoffs between number of enzymes, number of cisternae and enzyme
specificity to achieve given complexity

(a) Less complex target, 3-GMM approximation

(b) More complex target, 20-GMM approximation

(c) Less complex target, 3-GMM approximation

(d) More complex target, 20-GMM approximation

Figure 3: Tradeoffs amongst the glycan synthesis parameters, enzyme
specificity $\sigma$, cisternal number $N_{C}$ and enzyme number $N_{E}$, to
achieve a complex target distribution $\mathbf{c}^{\ast})$. (a)-(b) Normalised
Kullback-Leibler distance $\bar{D}(\sigma,N_{E},N_{C},\mathbf{c}^{\ast})$ as
function of $\sigma$ and $N_{C}$ (for fixed $N_{E}=3$), (c)-(d)
$\bar{D}(\sigma,N_{E},N_{C},\mathbf{c}^{\ast})$ as function of $\sigma$ and
$N_{E}$ (for fixed $N_{C}=3$), with the target distribution
$\mathbf{c}^{\ast}$ set to the $3$-GMM (less complex) and $20$-GMM (more
complex) approximations for the human T-cell MSMS data.
$\bar{D}(\sigma,N_{E},N_{C},\mathbf{c}^{\ast})$ is a convex function of
$\sigma$ for each $(N_{E},N_{C},\mathbf{c}^{\ast})$, decreasing in
$N_{C},N_{E}$ for each $(\sigma,\mathbf{c}^{\ast})$, increasing in the
complexity of $\mathbf{c}^{\ast}$ for fixed $(\sigma,N_{E},N_{C})$. The
specificity
$\sigma_{\min}(\mathbf{c}^{\ast},N_{E},N_{C})=\text{argmin}_{\sigma}\\{\bar{D}(\sigma,N_{E},N_{C},\mathbf{c}^{\ast})\\}$
that minimises the error for given $(N_{E},N_{C},\mathbf{c}^{\ast})$ is an
increasing function of $N_{C},N_{E}$ and the complexity of the target
distribution $\mathbf{c}^{\ast}$. Furthermore, the curvature of
$\bar{D}(\sigma,N_{E},N_{C},\mathbf{c}^{\ast})$ at
$\sigma_{\min}(N_{E},N_{C},\mathbf{c}^{\ast})$, related to sensitivity, is a
decreasing function of $N_{C},N_{E}$.

We are now in a position to catalogue the main results that follow from an
optimization of the parameters of the glycan synthesis machinery to a given
target distribution, Figs. 3-4

  1. 1.

The normalized KL-distance $\bar{D}(\sigma,N_{E},N_{C},\mathbf{c}^{\ast})$ is
a convex function of $\sigma$ for fixed values for other parameters (Fig. 3),
i.e. it first decreases with $\sigma$ and then increases beyond a critical
value of $\sigma_{\min}$. $\bar{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 reducing
the normalised distance $\bar{D}$ is approximately equal (Figs. 4a, 4b). The
lower complexity distributions can be synthesized with high fidelity with
small $(N_{E},N_{C})$, whereas higher complexity distributions require
significantly larger $(N_{E},N_{C})$, Figs. 4a, 4b. For a typical mammalian
cell, the number of enzymes in the N-glycosylation pathway are in the range
$N_{E}=10-20$ Umaña and Bailey (1997); Krambeck and Betenbaugh (2005);
Krambeck  _et al._ (2009); Fisher and Ungar (2016), Fig. 4b would then suggest
that the optimal cisternal number would range from $N_{C}=3-8$ Sengupta and
Linstedt (2011).

  2. 2.

The optimal enzyme specificity
$\sigma_{\min}(\mathbf{c}^{\ast},N_{C})=\text{argmin}_{\sigma}\\{\bar{D}(\sigma,\bar{N}_{E},N_{C},\mathbf{c}^{\ast})\\}$,
that minimises the error as function of $(N_{C},\mathbf{c}^{\ast})$ with
$N_{E}$ fixed at $\bar{N}_{E}$, is an increasing function of $N_{C}$ and the
complexity of the target distribution $\mathbf{c}^{\ast}$ (Figs. 3a, 3b, 4c,
4d). This is consistent with the results in Appendix E where we established
that the width of the synthesized distribution is inversely dependent on the
specificity $\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 $\bar{N}_{C}$, and $N_{E}$ is varied (Figs. 3c, 3d, 4c,
4d).

  3. 3.

Let $\sigma_{\min}(N_{C},N_{E},\mathbf{c}^{\ast})$ denote the value of
$\sigma$ that minimizes $\bar{D}(\sigma,N_{C},N_{E},\mathbf{c}^{\ast})$. Then
the second-derivative
$\nabla^{2}_{\sigma_{\min}}\bar{D}(N_{C},N_{E},\mathbf{c}^{\ast})=\frac{d^{2}}{d\sigma^{2}}\bar{D}(\sigma,N_{C},N_{E},\mathbf{c}^{\ast})\mid_{\sigma=\sigma_{\min}}$
denotes the curvature at $\sigma_{\min}$, and is measure of the sensitivity of
$\bar{D}(\sigma,N_{E},N_{C},\mathbf{c}^{\ast})$ to $\sigma$ for values close
to $\sigma_{\min}(N_{E},N_{C},\mathbf{c}^{\ast})$.
$\nabla^{2}_{\sigma_{\min}}\bar{D}(N_{C},N_{E},\mathbf{c}^{\ast})$ is a
decreasing function of $N_{C}$ (resp. $N_{E}$) for fixed values of
$(N_{E},\mathbf{c}^{\ast})$ (resp. $(N_{C},\mathbf{c}^{\ast})$), see Figs. 3,
4e, 4f. Thus, for any target distribution $\mathbf{c^{\ast}}$ there is a
minimal value of $(N_{E},N_{C})$ such that the target can be synthesized with
high fidelity provided the sensitivity $\sigma$ is tightly controlled at
$\sigma_{\min}(N_{C},N_{E},\mathbf{c}^{\ast})$, and there is larger value of
$(N_{E},N_{C})$ such that the target can be synthesized even if the control on
$\sigma$ is less tight.

Ungar et al. Fisher  _et al._ (2019) optimize incoming glycan ratio, transport
rate and effective reaction rates in order to synthesize a narrow target
distribution centred around a desired glycan. The ability to produce specific
glycans without much heterogeneity is an important goal in pharma industry.
They define heterogeneity as the total number of glycans synthesized, and show
that increasing the number of compartments $N_{C}$ decreases heterogeneity,
and increases the concentration of the specific glycan. They also show that
changing transport rate does not affect the heterogeneity. Our results are
entirely consistent with theirs - we have shown that $\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.

(a) Fidelity for less complex target, $\mathbf{c^{*}}=$ 3-GMM
approximation

(b) Fidelity for more complex target $\mathbf{c^{*}}=$ 20-GMM approximation

(c) Optimal enzyme specificity for less complex target, $\mathbf{c^{*}}=$
3-GMM approximation

(d) Optimal enzyme specificity for more complex target $\mathbf{c^{*}}=$
20-GMM approximation

(e) Sensitivity for less complex target, $\mathbf{c^{*}}=$ 3-GMM approximation

(f) Sensitivity for more complex target $\mathbf{c^{*}}=$ 20-GMM approximation

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

###  VI.3 Optimal partitioning of enzymes in cisternae

(a)

(b)

Figure 5: Optimal enzyme partitioning in cisternae. (a) Heat map of the
effective reaction rates in each cisterna (representing the optimal enzyme
partitioning) and the steady state concentration in the last compartment
($\mathbf{c}^{(N_{C})}$) for the 20-GMM target distribution. Here $N_{E}=5$,
$N_{C}=7$, normalised $D_{KL}({\bf T}^{(20)}\|{\bf c}^{(N_{C})})/H({\bf
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.

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 our
chemical reaction model (Sect. IV.1) requires some care, since it incorporates
the following enzymatic features: (a) enzymes with a finite specificity
$\sigma$ can catalyse several reactions, although with an efficiency that
varies with both the substrate index $k$ and cisternal index $j$, and (b)
every enzyme appears in each cisternae; however their reaction efficiencies
depend on the enzyme levels, the enzymatic reaction rates and the enzyme
matching function $\bf{L}$, all of which depend on the cisternal index $j$.

Thus, rather than determining the cisternal partitioning of enzymes, we
instead identify chemical reactions that occur with high propensity in each
cisternae. For this we define an effective reaction rate $\bar{R}(j,k)$ for
${\cal{P}}c_{k}\to{\cal{P}}c_{k+1}$ in the $j$-th cisterna as

${\bar{R}}(j,k)=\sum_{\alpha=1}^{N_{E}}R^{(j)}_{\alpha}P^{(j)}(k,\alpha).$
(10)

According to our model presented in Sect. IV.1, the list of reactions with
high effective reaction rates in each cisterna, corresponds to a cisternal
partitioning of the perfect enzymes. In a future study, we will consider a
Boolean version of a more complex chemical model, to address more clearly, the
optimal enzyme partitioning amongst cisternae.

Figure 5 (a) (i) shows the heat map of the effective reaction rates in each
cisterna for the optimal $N_{E},N_{C},\sigma$ that minimises the normalised
KL-distance to the 20 GMM target distribution $\mathbf{T}^{(20)}$ (see Fig. 5
(a) (ii)). The optimized glycan profile displayed in Fig. 5 (a) (iii) is very
close to the target. An interesting observation from Fig. 5a(i) is that the
same reaction can occur in multiple cisternae.

Keeping everything else fixed at the optimal value, we ask whether simply
repartitioning the optimal enzymes amongst the cisternae, alters the displayed
glycan distribution. In Fig. 5 (b) (i), we have exchanged the enzymes of the
fourth and second cisterna. The glycan profile after enzyme partitioning (see
Fig. 5 (b) ((iii)) is now completely altered (compare Fig. 5 (b) (ii) with
Fig. 5 (b) (iii)). Thus one may achieve a different glycan distribution by
repartitioning enzymes amongst the same number of cisternae Jaiman and Thattai
(2018).

##  VII Strategies to achieve high glycan diversity

(a)

(b)

Figure 6: Strategies for achieving high glycan diversity. Diversity versus
$N_{C}$ and transport rate $\mu$ at various values of specificity $\sigma$ for
fixed $N_{E}=3$. (a) Diversity vs. $N_{C}$ at optimal transport rate $\mu$.
Diversity initially increases with $N_{C}$, but eventually levels off. The
levelling off starts at a higher $N_{C}$ when $\sigma$ is increased. These
curves are bounded by the $\sigma=0$ curve. (b) Diversity vs. cisternal
residence time ($\mu^{-1}$) in units of the reaction time ($R_{\min}^{-1}$) at
various value of $\sigma$, for fixed $N_{C}=4$ and $N_{E}=10$. This has
implications for glycoengineering (see text) where the task is to produce a
particular glycan profile with low heterogeneity Fisher  _et al._ (2019);
Jaiman and Thattai (2018).

So far we have studied how the complexity of the target glycan distribution
places constraints on the evolution of Golgi cisternal number and enzyme
specificity. We now take up another issue, namely, how the physical properties
of the Golgi cisternae, namely cisternal number and inter-cisternal transport
rate, may drive diversification of glycans 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).

In this section, we study how changing the physical parameters in our chemical
synthesis model can lead to changes in the diversity of glycan distributions.

We define diversity as the total number of glycan species produced above a
specified threshold abundance $c_{th}$. This last condition is necessary
because very small peaks will not be distinguishable in the presence of noise.
In computing the diversity from our chemical synthesis model, we have chosen
the threshold to be $c_{th}=1/N_{s}$, where $N_{s}$ is the total number of
glycan species. We have checked that the qualitative results do not depend on
this choice, Fig. A6.

Using the sigmoid function $(1+e^{-x/\tau})^{-1}$ as a continuous
approximation to the Heaviside function $\Theta(x)$, we define the following
optimization to achieve the maximal diversity for a given set of parameter
values, $N_{E},N_{C},\sigma$,

$\displaystyle\text{Diversity}(\sigma,N_{C},N_{E}):=$max_μ, R, L∑_i=1^N_s (1+
e^-N_s(c_i - c_th))^-1s.t.R_min ≤R^(j)_α≤R_max,μ_min ≤μ^(j) ≤μ_max, (13)

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 F for details on the parameter estimation.

The results displayed in Fig. 6(a), 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 correlative at best, since there are
many ways to achieve high glycan diversity - notably by increasing the number
of enzymes.

On the other hand, one of the goals of glycoengineering is to produce a
particular glycan profile with low heterogeneity 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 early exit of incomplete oligomers. At
large cisternal residence time the diversity again decreases as more reactions
are taken to completion. Note that the peak is generally very flat, this is
consistent with the results of 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.

##  VIII Discussion

The precision of the stereochemistry and enzymatic kinetics of these
N-glycosylation reactions Varki  _et al._ (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 may help predict
glycosylation patterns and direct glycoengineering projects to optimize
glycoform distributions.

In this paper, we have been interested in the role of glycans as a marker or
molecular code of cell identity Gabius (2018); Varki (2017); Pothukuchi  _et
al._ (2019). In particular, we have studied one aspect of molecular coding,
namely the fidelity of this glycan code generated by enzymatic and transport
processes located in the secretory apparatus of the cell. This involves a
method of analysis that draws on many different fields, and so it might be
useful to provide a short summary of the assumptions, methods and results of
the paper:

  1. 1.

The distribution of glycans at the cell surface is a marker of cell-type
identity Varki  _et al._ (2009); Gabius (2018); Varki (2017); Pothukuchi  _et
al._ (2019). This glycan distribution can be very complex; it is believed that
there is an evolutionary drive for having glycan distributions of high
complexity arising from the following considerations,

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

    2. (b)

Consequence of pathogen-mediated selection pressures Varki  _et al._ (2009);
Varki (2017); Gagneux and Varki (1999).

    3. (c)

Consequence of herd immunity within a heterogenous population of cells of a
community WILLS and GREEN (1995) or within a single organism Drickamer and
Taylor (1998).

  2. 2.

The glycans at the cell surface are the end product of a sequential chemical
processing via a set of enzymes resident in the Golgi cisternae, and transport
across cisternae Varki (2017, 1998); Pothukuchi  _et al._ (2019). Using a
fairly general and tractable model for chemical synthesis and transport, we
compute the synthesized glycan distribution at the cell surface. Parameters of
our synthesis model include the number of enzymes $N_{E}$, specificity of
enzymes $\sigma$, number of cisternae $N_{C}$ and transport rates $\mu$.

  3. 3.

We measure the reliability or fidelity of cell identity Stanley (2011);
Pothukuchi  _et al._ (2019); Varki (1998) in terms of the error between
synthesized glycan distribution on the cell surface from the its internal
“target” distribution using the Kullback- Leibler distance $D_{KL}$ Cover and
Thomas (2012); MacKay (2003). In our numerical study, we obtain the target
distribution for the given cell type by suitable coarse-graining of the MSMS
data for the human T-cells Cummings and Crocker (2020). We solve a constrained
optimization problem for minimising $D_{KL}$, and study the tradeoffs between
$N_{E}$, $N_{C}$ and $\sigma$.

  4. 4.

The results of the optimization to achieve a given target complexity are
summarised in Figs. 3-4. Here, we highlight some its direct consequences:

    1. (a)

Keeping the number of enzymes fixed, a more elaborate transport mechanism (via
control of $N_{C}$ and $\mu$) is essential for synthesising high complexity
target distributions (Figs. 4a, 4b). Fewer cisternae cannot be compensated for
by optimising the enzymatic synthesis (via control of parameters
$\mathbf{R,L}$ and $\sigma$).

    2. (b)

Thus, our study suggests that fidelity of the glycan code generation provides
a functional control of Golgi cisternal number. It also provides a
quantitative argument for the evolutionary requirement of multiple-
compartments, by demonstrating that the fidelity and sensitivity of the glycan
code arising from a chemical synthesis that involves multiple cisternae is
higher than one that involves a single cisterna (keeping everything else
fixed) (Figs. 4a, 4b, 4e, 4f). This feature that with multiple cisternae and
precise enzyme partitioning, one may generically achieve a highly accurate
representation of the target distribution, has been highlighted in an
algorithmic model of glycan synthesis Jaiman and Thattai (2018).

    3. (c)

Our study shows that for a fixed $N_{C}$ and $N_{E}$, there is an optimal
enzyme specificity that achieves the lowest distance from a given target
distribution. As we see in Fig. 4d, this optimal enzyme specificity can be
very high for highly complex target distributions.

    4. (d)

Organisms such as plants and algae, have a diverse repertoire of glycans that
are utilised 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 (Fig. 6). However, this compromises on the complexity of the glycan
distribution, revealing a tension between complexity and diversity. One way of
relieving this tension is to have larger $N_{E}$ and $N_{C}$.

    5. (e)

Consider a situation where the environment, and hence the target glycan
distribution, fluctuates rapidly. When synthesis parameters cannot change
rapidly enough to track the environment, high specificity enzymes can lead to
a lowering of the cell’s fitness Nam  _et al._ (2012); Peracchi (2018). Having
slightly sloppy enzymes may give the best selective advantage in a time
varying environment. This compromise, between robustness in a changing
environment and the demand for complexity, is achieved by having sloppy
enzymes, that allows the system to be more evolvable Nam  _et al._ (2012);
Peracchi (2018). However, sloppy enzymes are subject to errors from
synthesising the wrong reaction products. In this case, error correcting
mechanisms must be in place to ensure fidelity of the glycan code. We leave
the role of intra-cellular transport in providing non-equilibrium proof-
reading mechanisms to reduce such coding errors, and its optimal adaptive
strategies and plasticity in a time varying environment, as a task for the
future.

Admittedly the chemical network that we have considered here is much simpler
than the chemical network associated with all possible protein modifications
in the secretory pathway. For instance, typical N-glycosylation pathways would
involve the glycosylation of a variety of GBPs. Further, apart from
N-glycosylation, there are other glycoprotein, proteoglycan and glycolipid
synthesis pathways Alberts  _et al._ (2002); Varki  _et al._ (2009);
Pothukuchi  _et al._ (2019). We believe our analysis is generalisable and that
the qualitative results we have arrived at would still hold.

To conclude, our work establishes the link between the cisternal machinery
(chemical and transport) and optimal coding. We find that the pressure to
achieve the target glycan code for a given cell type, places strong
constraints on the cisternal number and enzyme specificity Sengupta and
Linstedt (2011). An important implication is that a description of the
nonequilibrium self-assembly of a fixed number of Golgi cisternae must combine
the dynamics of chemical processing and membrane dynamics involving fission,
fusion and transport Sengupta and Linstedt (2011); Sachdeva  _et al._ (2016);
Sens and Rao (2013). We believe this is a promising direction for future
research.

##  IX Acknowledgments

We thank M. Thattai, A. Jaiman, S. Ramaswamy, 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 for consultations on the MSMS data and literature. We
acknowledge the computational facilities at NCBS. MR acknowledges a JC Bose
Fellowship from DST (Government of India), and thanks Institut Curie for
hosting a visit under the Labex program. 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. QV thanks the Simons Centre (NCBS) for hosting his
visit.

Appendix

##  Appendix A Kinetics of sequential chemical reactions and transport

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 substrate concentrations $c_{k}^{(j)}$ change with time as,

$\begin{split}\frac{dc^{(1)}_{1}}{dt}=\;&q-\sum_{\alpha=1}^{N_{E}}\frac{V(1,1,\alpha)P^{(1)}(1,\alpha)c_{1}^{(1)}}{M(1,1,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(1)}(k^{\prime},\alpha)c_{k^{\prime}}^{(1)}}{M(1,k^{\prime},\alpha)}\right)}-\mu^{(1)}c_{1}^{(1)}\\\
\frac{dc_{k}^{(1)}}{dt}=&\sum_{\alpha=1}^{N_{E}}\frac{V(1,k-1,\alpha)P^{(1)}(k-1,\alpha)c_{k-1}^{(1)}}{M(1,k-1,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(1)}(k^{\prime},\alpha)c_{k^{\prime}}^{(1)}}{M(1,k^{\prime},\alpha)}\right)}\\\
&\mbox{}-\sum_{\alpha=1}^{N_{E}}\frac{V(1,k,\alpha)P^{(1)}(k,\alpha)c_{k}^{(1)}}{M(1,k,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(1)}(k^{\prime},\alpha)c_{k^{\prime}}^{(1)}}{M(1,k^{\prime},\alpha)}\right)}\mbox{}-\mu^{(1)}c_{k}^{(1)}\\\
\frac{dc^{(1)}_{N_{s}}}{dt}=&\sum_{\alpha=1}^{N_{E}}\frac{V(1,N_{s}-1,\alpha)P^{(1)}(N_{s}-1,\alpha)c_{N_{s}-1}^{(1)}}{M(1,N_{s}-1,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(1)}(k^{\prime},\alpha)c_{k^{\prime}}^{(1)}}{M(1,k^{\prime},\alpha)}\right)}-\mu^{(1)}c_{N_{s}}^{(1)}\\\
\end{split}$ (14)

for cisterna-$1$, and

$\begin{split}\frac{dc^{(j)}_{1}}{dt}=\;&\mu^{(j-1)}c_{1}^{(j-1)}-\sum_{\alpha=1}^{N_{E}}\frac{V(j,1,\alpha)P^{(j)}(1,\alpha)c_{1}^{(j)}}{M(j,1,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(j)}(k^{\prime},\alpha)c_{k^{\prime}}^{(j)}}{M(j,k^{\prime},\alpha)}\right)}-\mu^{(j)}c^{(j)}_{1}\\\
\frac{dc_{k}^{(j)}}{dt}=\;&\mu^{(j-1)}c_{k}^{(j-1)}+\mbox{
}\sum_{\alpha=1}^{N_{E}}\frac{V(j,k-1,\alpha)P^{(j)}(k-1,\alpha)c_{k-1}^{(j)}}{M(j,k-1,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(j)}(k^{\prime},\alpha)c_{k^{\prime}}^{(j)}}{M(j,k^{\prime},\alpha)}\right)}\\\
&\mbox{}-\sum_{\alpha=1}^{N_{E}}\frac{V(j,k,\alpha)P^{(j)}(k,\alpha)c_{k}^{(j)}}{M(j,k,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(j)}(k^{\prime},\alpha)c_{k^{\prime}}^{(j)}}{M(j,k^{\prime},\alpha)}\right)}-\mu^{(j)}c_{k}^{(j)}\\\
\frac{dc_{N_{s}}^{(j)}}{dt}=\;&\mu^{(j-1)}c_{N_{s}}^{(j-1)}+\sum_{\alpha=1}^{N_{E}}\frac{V(j,N_{s}-1,\alpha)P^{(j)}(N_{s}-1,\alpha)c_{N_{s}-1}^{(j)}}{M(j,N_{s}-1,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(j)}(k^{\prime},\alpha)c_{k^{\prime}}^{(j)}}{M(j,k^{\prime},\alpha)}\right)}-\mu^{(j)}c_{N_{s}}^{(j)}\end{split}$
(15)

for cisternae $j=2,3,\ldots,N_{C}$. These set of dynamical equations
(14)-(15), with initial conditions, can be solved to obtain the concentration
$c^{(j)}_{k}(t)$ for $t\geq 0$.

Equations (14)-(15) automatically obey the conservation law for the protein
concentration ($p$), i.e., denoting the protein concentration in the $j$-th
cisterna as $p^{(j)}=\sum_{k^{\prime}=1}^{N_{s}}c_{k^{\prime}}^{(j)}$, we
automatically obtain,

$\displaystyle\frac{dp^{(1)}}{dt}$ $\displaystyle=$ $\displaystyle
q-\mu^{(1)}p^{(1)}$ $\displaystyle\frac{dp^{(j)}}{dt}$ $\displaystyle=$
$\displaystyle\mu^{(j-1)}p^{(j-1)}-\mu^{(j)}p^{(j)}$

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

At steady state, the left hand side of the above equations is set to zero,
which after rescaling, gives the nonlinear recursion relations displayed in
(5) and (6) of the main text.

##  Appendix B A computationally simpler optimization equivalent to
Optimization A

Define a new set of parameters,

$R(j,k,\alpha)=\sum_{\alpha=1}^{N_{E}}\frac{V(j,k,\alpha)}{M(j,k,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(j)}(k^{\prime},\alpha)c_{k^{\prime}}^{(j)}}{M(j,k^{\prime},\alpha)}\right)}$
(16)

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

$\begin{split}v_{1}^{(1)}=&\frac{1}{\mu^{(1)}+\sum_{\alpha=1}^{N_{E}}R(1,1,\alpha)P^{(1)}(1,\alpha)}\\\
v_{k}^{(1)}=&\frac{v_{k-1}^{(1)}\sum_{\alpha=1}^{N_{E}}R(1,k-1,\alpha)P^{(1)}(k-1,\alpha)}{\mu^{(1)}+\sum_{\alpha=1}^{N_{E}}R(1,k,\alpha)P^{(1)}(k,\alpha)}\\\
v_{N_{s}}^{(1)}=&\frac{v_{N_{s}-1}^{(1)}\sum_{\alpha=1}^{N_{E}}R(1,N_{s}-1,\alpha)P^{(1)}(N_{s}-1,\alpha)}{\mu^{(1)}}\end{split}$
(17)

for $j=1$, and

$\begin{split}v_{1}^{(j)}=\;&\frac{v_{1}^{(j-1)}\mu^{(j-1)}}{\mu^{(j)}+\sum_{\alpha=1}^{N_{E}}R(j,1,\alpha)P^{(j)}(1,\alpha)}\\\
v_{k}^{(j)}=\;&\frac{v_{k}^{(j-1)}\mu^{(j-1)}}{\mu^{(j)}+\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)P^{(j)}(k,\alpha)}\\\
&\mbox{}+\frac{v_{k-1}^{(j)}\sum_{\alpha=1}^{N_{E}}R(j,k-1,\alpha)P^{(j)}(k-1,\alpha)}{\mu^{(j)}+\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)P^{(j)}(k,\alpha)}\\\
v_{N_{s}}^{(j)}=\;&\frac{v_{N_{s}}^{(j-1)}\sum_{\alpha=1}^{N_{E}}R(j,N_{s}-1,\alpha)P^{(j)}(N_{s}-1,\alpha)}{\mu^{(j)}}+\frac{v_{N_{s}}^{(j-1)}\mu^{(j-1)}}{\mu^{(j)}}\end{split}$
(18)

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

In Appendix C we show that for $\mathbf{v}$ obtained from (17)-(18) for any
parameter $(\mathbf{R},\mathbf{L})$, there exists parameter $(\mathbf{M,V,L})$
such that (5)-(6) are automatically satisfied when we set
$\mathbf{c}=\mathbf{v}$, i.e. $\mathbf{v}$ is the steady state concentration
for $(\mathbf{M},\mathbf{V},\mathbf{L})$. Thus, the set of all concentration
profiles defined by (17)-(18) as a function of all possible values of the
parameters $(\mathbf{R},\mathbf{L})$ is identical to the set defined by
(5)-(6) as function of $(\mathbf{M},\mathbf{V},\mathbf{L})$. This is a crucial
insight, since it allows us to search the entire parameter space using
(17)-(18), where the concentration is known explicitly in terms of
$(\mathbf{R},\mathbf{L})$. See Figure A1 for a flow chart of the two
optimization schemes.

To pose this new optimization problem, it is convenient to define
$\bar{v}_{i}=\mu^{(N_{c})}v^{(N_{c})}_{i}$. Then, it follows that

Optimization B:

$D(\sigma,N_{E},N_{C},\mathbf{c}^{\ast}):=\min_{\mathbf{R}\geq\mathbf{0},\
\mathbf{L}}D_{KL}(\mathbf{c}^{\ast}\|\bar{\mathbf{v}})$ (19)

is equivalent to (8). Since $\mathbf{v}$ is explicitly known as a function of
$(\mathbf{R},\mathbf{L})$, optimization B (19) is a more tractable
optimization problem than (8). Note that in this setting, the function
$D(\sigma,N_{E},N_{C},\mathbf{c}^{\ast})$ (19) is independent of the rates
$\bm{\mu}$.

While this optimization is easy to implement, we note that the parameters
(e.g., reaction rates, specificity) are not constrained to take only
physically relevant values; a legitimate concern is that the absence of such
physico-chemical constraints might drive this optimization to physically
unrealistic solutions.

There are two possible ways to impose these parameter constraints. One is to
impose constraints on the “microscopic” chemical parameters, such as the rate
of individual reactions $R(j,k,\alpha)$ and the inter-cisternal transport rate
${\mu}^{(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.

Optimization C :

$\displaystyle D(\sigma,N_{C},N_{E},{\bf{c^{*}}}):=$min_μ, R, LD_KL(c^∗∥
¯v)s.t.R_min ≤R(j,k,α) ≤R_max,μ_min ≤μ^(j) ≤μ_max. (22)

The upper and lower bounds on the rates $\mathbf{R}$ and $\bm{\mu}$ are
estimated in Appendix F : $\mu_{\max}=1$/min (resp. $\mu_{\min}=.01$/min) and
$R_{\max}=20$/min (resp. $R_{\min}=.018$/min).

Figure A1: 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
$\bf{c}$ is equal to the set of all $\bf{v}$. We additionally establish that
the optimum $\mathbf{v}_{\min}=\mathbf{c}_{\min}$.

##  Appendix C Equivalence of Optimizations $A$ and $B$

Let

$\mathcal{A}=\left\\{[c_{k}^{(j)}]_{j,k}:\begin{array}[]{l}M(j,k,\alpha)\geq
0,V(j,k,\alpha)\geq 0,l_{\alpha}^{(j)}\geq 0,\\\ \text{$[c_{k}^{(j)}]_{jk}$
given by \eqref{eq:concentrations1} and
\eqref{eq:concentrationsj}}\end{array}\right\\}$

denote the set of concentrations achievable in Optimization $A$. Similarly,
let

$\mathcal{B}=\left\\{[v_{k}^{(j)}]_{j,k}:\begin{array}[]{l}R(j,k,\alpha)\geq
0,l_{\alpha}^{(j)}\geq 0\\\ \text{$[v_{k}^{(j)}]_{j,k}$ given by
\eqref{eq:Bconcentrations1} and
\eqref{eq:Bconcentrationsj}}\end{array}\right\\}$

denote the set of concentrations achievable in the Optimization $B$.

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

$R(j,k,\alpha)=\sum_{\alpha=1}^{N_{E}}\frac{V(j,k,\alpha)}{M(j,k,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(j)}(k^{\prime},\alpha)c_{k^{\prime}}^{(j)}}{M(j,k^{\prime},\alpha)}\right)}\;\geq
0$

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

Suppose $[v_{k}^{(j)}]_{j,k}\in\mathcal{B}$. Let $[R(j,k,\alpha)]$,
$[l_{\alpha}^{(j)}]$ denote the corresponding parameters. Since
$\sum_{k=1}^{N_{s}}v_{k}^{(j)}=1/\mu^{(j)}<\infty$, it follows that
$\sum_{k=1}^{N_{s}}P^{(j)}(k,\alpha)v_{k}^{(j)}<1/\mu^{(j)}<\infty$. Thus,
there exists parameters $[M(j,k,\alpha)]$, $[V(j,k,\alpha)]$ and
$[l_{\alpha}^{(j)}]$ such that

$R(j,k,\alpha)=\sum_{\alpha=1}^{N_{E}}\frac{V(j,k,\alpha)}{M(j,k,\alpha)\left(1+\sum_{k^{\prime}=1}^{N_{s}}\frac{P^{(j)}(k^{\prime},\alpha)v_{k^{\prime}}^{(j)}}{M(j,k^{\prime},\alpha)}\right)}$
(23)

Therefore, $[v_{k}^{(j)}]_{j,k}$ satisfy (5) and (6), i.e.
$[v_{k}^{(j)}]_{j,k}\in\mathcal{A}$.

Moreover, suppose $\mathbf{v}$ satisfies (17)-(18) for a given set of
parameters $(\mathbf{R},\mathbf{L})$. Then there exist
$(\mathbf{M},\mathbf{V},\mathbf{L})$ such that $\mathbf{v}$ satisfies (5)-(6),
i.e. $\mathbf{v}$ is the steady state concentration for
$(\mathbf{M},\mathbf{V},\mathbf{L})$.

##  Appendix D Analytical solution when $N_{s}\gg 1$

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

$\displaystyle\frac{dc^{(1)}_{k}}{dk}\approx c^{(1)}_{k}-c^{(1)}_{k-1}$
$\displaystyle=\left(\frac{\sum_{\alpha=1}^{N_{E}}R(1,k-1,\alpha)\exp(-\sigma|k-1-l_{\alpha}^{(1)}|)}{\mu^{(1)}+\sum_{\alpha=1}^{N_{E}}R(1,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(1)}|)}-1\right)c^{(1)}_{k-1}$
$\displaystyle\approx-\left(\frac{\mu^{(1)}+\frac{d}{dk}\sum_{\alpha=1}^{N_{E}}R(1,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(1)}|)}{\mu^{(1)}+\sum_{\alpha=1}^{N_{E}}R(1,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(1)}|)}\right)c^{(1)}_{k},$

and

$\displaystyle\frac{dc^{(j)}_{k}}{dk}\approx c^{(j)}_{k}-c^{(j)}_{k-1}$
$\displaystyle=$
$\displaystyle\frac{\mu^{(j-1)}}{\mu^{(j)}+\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(j)}|)}c^{(j-1)}_{k}$
$\displaystyle\mbox{}-\left(\frac{\mu^{(j)}+\frac{d}{dk}\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(j)}|)}{\mu^{(j)}+\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(j)}|)}\right)c^{(j)}_{k}$

for $j=2,\ldots,N_{C}$. In (LABEL:DifferentialEquantion1) and
(LABEL:DifferentialEquationNc),

$\displaystyle\frac{d}{dk}\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(j)}|)$
(26)
$\displaystyle=\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)\sigma\exp(-\sigma|k-l_{\alpha}^{(j)}|)(1-2\mathbb{I}(k\geq
l_{\alpha}))+\;R^{\prime}(j,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(j)}|)$

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

Define a vector function $C(k)\in\mathbb{R}^{N}_{c}$ of the continuous
variable $k$ by $C(k)=[c^{(1)}_{k},c^{(2)}_{k},\ldots c^{(N_{C})}_{k}]$. Then
(LABEL:DifferentialEquantion1) and (LABEL:DifferentialEquationNc) can be
written as:

$\frac{dC(k)}{dk}=M(k)C(k)$ (27)

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

$\displaystyle M(k)$ $\displaystyle=$
$\displaystyle\begin{bmatrix}A^{(1)}(k)&0&0&0&\ldots 0\\\
B^{(2)}(k)&A^{(2)}(k)&0&0&\ldots 0\\\ 0&B^{(3)}(k)&A^{(3)}(k)&0&\ldots 0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&\ldots&0&B^{({N_{C}})}(k)&A^{({N_{C}})}(k)\\\ \end{bmatrix}$

with

$\displaystyle A^{(j)}(k)$ $\displaystyle=$
$\displaystyle-\frac{\mu^{(j)}+\frac{d}{dk}\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(j)}|)}{\mu^{(j)}+\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(j)}|)}$
$\displaystyle B^{(j)}(k)$ $\displaystyle=$
$\displaystyle\frac{\mu^{(j-1)}}{\mu^{(j)}+\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(j)}|)}$

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

$C(k)=C_{0}\exp\left(\Omega(k)\right)$ (29)

is written in terms of the Magnus Function
$\Omega(k)=\sum_{n=1}^{\infty}\Omega(n,k)$, obtained from the Baker-Campbell-
Hausdorff formula Blanes  _et al._ (2009),

$\displaystyle\Omega(1,k)$ $\displaystyle=$
$\displaystyle\int_{0}^{k}M(k_{1})dk_{1}$ $\displaystyle\Omega(2,k)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\int_{0}^{k}dk_{1}\int_{0}^{k_{1}}dk_{2}\left[M(k_{1}),M(k_{2})\right]$
$\displaystyle\Omega(3,k)$ $\displaystyle=$
$\displaystyle\frac{1}{6}\int_{0}^{k}dk_{1}\int_{0}^{k_{1}}dk_{2}\int_{0}^{k_{2}}dk_{3}\left[M(k_{1}),\left[M(k_{2}),M(k_{3})\right]\right]\mbox{}+\left[M(k_{3}),\left[M(k_{2}),M(k_{1})\right]\right]$
$\displaystyle\ldots\ldots$

where $\left[M(k_{1}),M(k_{2})\right]:=M(k_{1})M(k_{2})-M(k_{2})M(k_{1})$ is
the commutator, and the higher order terms in $\ldots$ contain higher order
nested commutators.

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

The commutator

$\displaystyle[M(k_{1}),M(k_{2})]=$[0 0 0 0 …0 a210 0 0 …0 a31a320 0 …0 0
a42a430 …0 ⋮⋮⋮⋮⋮0 …an,n-2an,n-10]

where

$\displaystyle a_{i,i-1}$ $\displaystyle=$ $\displaystyle
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})$
$\displaystyle a_{i,i-2}$ $\displaystyle=$ $\displaystyle
B^{(i-1)}(k_{2})B^{(i)}(k_{1})-B^{(i-1)}(k_{1})B^{(i)}(k_{2})$

The general form of $\Omega(n,k)$ is given by Blanes  _et al._ (2009)

$\begin{split}\Omega(n,k)=&\frac{z_{n}}{n!}\int_{0}^{k}dk_{1}\int_{0}^{k_{1}}dk_{2}\ldots\int_{0}^{k_{n-2}}dk_{n-1}\int_{0}^{k_{n-1}}dk_{n}\sum_{l}W_{l}M(k_{p^{l}_{1}})M(k_{p^{l}_{2}})\ldots
M(k_{p^{l}_{n}})\end{split}$ (30)

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

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

$\bar{M}=\begin{bmatrix}\bar{A}&0&0&0\ldots&0\\\
\bar{A}&\bar{A}&0&0\ldots&0\\\ 0&\bar{A}&\bar{A}&0\ldots&0\\\
\vdots&\vdots&\vdots&&\vdots\\\ 0&\ldots&0&\bar{A}&\bar{A}\\\ \end{bmatrix}$

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

$\displaystyle\Omega_{lm}(n,k)$ $\displaystyle\leq$ $\displaystyle
z_{n}\bar{M}^{n}_{l,m}\int_{0}^{k}dk_{1}\int_{0}^{k_{1}}dk_{2}\ldots\int_{0}^{k_{n-1}}dk_{n}$
(31) $\displaystyle=$ $\displaystyle
z_{n}\bar{M}^{n}\Big{|}_{lm}\frac{k^{n}}{n!}$

The matrix

$\bar{M}^{n}=\begin{bmatrix}a_{11}&0&0&0&\ldots&0\\\
a_{21}&a_{22}&0&0&\ldots&0\\\ a_{31}&a_{32}&a_{33}&0&\ldots&0\\\
a_{41}&a_{42}&a_{43}&a_{44}&\ldots&0\\\ \vdots&\vdots&\vdots&\vdots&&\vdots\\\
a_{n1}&\ldots&&a_{n,n-2}&a_{n,n-1}&a_{nn}\\\ \end{bmatrix}$

where $a_{lm}=S_{lm}(n)\bar{A}^{n}$ for appropriately defined polynomials
$S_{l,m}(n)$. Thus, it follows that $\Omega_{lm}\leq
z_{n}S_{lm}(n)(A^{*})^{n}\frac{k^{n}}{n!}$ and
$\Omega_{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$, i.e.
$k\leq\frac{1}{\bar{A}}$. Assuming $\mu^{(j)}=\mu\;\forall j$, we can bound
$\bar{A}$ as

$\displaystyle\bar{A}\leq\max_{j,k}\left(\frac{\mu+\sigma\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(j)}|)}{\mu+\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(j)}|)}\right.+\left.\frac{\sum_{\alpha=1}^{N_{E}}R^{\prime}(j,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(j)}|)}{\mu+\sum_{\alpha=1}^{N_{E}}R(j,k,\alpha)\exp(-\sigma|k-l_{\alpha}^{(j)}|)}\right)$
(32)

Since the parameters $\mu,\;\sigma,\;R(j,k,\alpha),\;l^{(j)}_{\alpha}$ and
$N_{E}$ are finite and positive, and $R^{\prime}(j,k,\alpha)$ is finite,
$\bar{A}$ has a finite upper bound, implying $k$ is always greater than zero,
and the series has finite radius of convergence.

While in principle we can obtain the glycan profile for any $N_{E}$ and
$N_{C}$ with arbitrary accuracy, assuming $R(j,k,\alpha)=R_{\alpha}^{(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)=\left\\{\begin{array}[]{ll}c_{0}e^{-k}\left(\frac{\mu+R\exp(-\sigma(l-k))}{\mu+R\exp(-\sigma
l)}\right)^{(1/\sigma)-1}&k\leq l\\\
c(l)e^{-(k-l)}\left(\frac{\mu+R}{\mu+R\exp(-\sigma(k-l))}\right)^{(1/\sigma)+1}&k>l\end{array}\right.$
(33)

A representative concentration profile is plotted in Fig. A2(a). The
concentration profile consists of two distinct components: an initial
exponential decay, and then an exponential rise and fall concentrated around
$l$. The relative weight of these two components is controlled by the
sensitivity $\sigma$ and the rate $R$. Such explicit formulae can be obtained
for any $N_{E}>1$, as long as $N_{C}=1$.

(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{split}c^{(2)}(k)=\;&c_{0}\mu^{(1)}e^{-k}\Big{(}\frac{\mu^{(2)}+R^{(2)}\exp(-\sigma(l^{(2)}-k))}{\mu^{(1)}+R^{(1)}e^{-\sigma
l^{(1)}}}\Big{)}^{(1/{\sigma})-1}\int_{0}^{k}\frac{(\mu^{(1)}+R^{(1)}\exp(-\sigma(l^{(1)}-k)))^{(1/\sigma)-1}}{(\mu^{(2)}+R^{(2)}\exp(-\sigma(l^{(2)}-k)))^{1/\sigma}}dk\\\
&\mbox{}+c^{(2)}(0)e^{-k}\left(\frac{\mu^{(2)}+R^{(2)}e^{-\sigma(l^{(2)}-k)}}{\mu^{(2)}+R^{(2)}e^{-\sigma
l^{(2)}}}\right)^{(1/\sigma)-1}\end{split}$ (34)

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

$\displaystyle c^{(2)}(k)$ $\displaystyle=$ $\displaystyle
c^{(1)}(l^{(1)})\mu^{(1)}e^{-(k-l^{(1)})}(\mu^{(1)}+R^{(1)})^{(1/\sigma)+1}(\mu^{(2)}+R^{(2)}\exp(-\sigma(l^{(2)}-k)))^{(1/\sigma)-1}$
(35)
$\displaystyle\int_{l^{(1)}}^{k}\frac{(\mu^{(2)}+R^{(2)}\exp(-\sigma(l^{(2)}-k)))^{-1/\sigma}}{(\mu^{(1)}+R^{(1)}\exp(-\sigma(k-l^{(1)})))^{(1/\sigma)+1}}dk$
$\displaystyle+$ $\displaystyle
c^{(2)}(l^{(1)})e^{-(k-l^{(1)})}\left(\frac{\mu^{(2)}+R^{(2)}e^{-\sigma(l^{(2)}-k)}}{\mu^{(2)}+R^{(2)}e^{-\sigma(l^{(2)}-l^{(1)})}}\right)^{(1/\sigma)-1}$

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

$\displaystyle c^{(2)}(k)$ $\displaystyle=$ $\displaystyle
c^{(1)}(l^{(1)})\mu^{(1)}e^{-(k-l^{(1)})}\left(\frac{\mu^{(1)}+R^{(1)}}{\mu^{(2)}+R^{(2)}\exp(-\sigma(k-l^{(2)}))}\right)^{(1/\sigma)+1}$
(36)
$\displaystyle\int_{l^{(2)}}^{k}\frac{(\mu^{(2)}+R^{(2)}\exp(-\sigma(k-l^{(2)})))^{1/\sigma}}{(\mu^{(1)}+R^{(1)}\exp(-\sigma(k-l^{(1)})))^{(1/\sigma)+1}}dk$
$\displaystyle+$ $\displaystyle
c^{(2)}(l^{(2)})e^{-(k-l^{(2)})}\left(\frac{\mu^{(2)}+R^{(2)}}{\mu^{(2)}+R^{(2)}e^{-\sigma(k-l^{(2)})}}\right)^{(1/\sigma)+1}$

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

$\displaystyle c^{(2)}(k)$ $\displaystyle=$ $\displaystyle
c_{0}\mu^{(1)}e^{-k}\frac{(\mu^{(1)}+R^{(1)}e^{-\sigma
l^{(1)}})^{1-(1/\sigma)}}{(\mu^{(2)}+R^{(2)}\exp(-\sigma(k-l^{(2)})))^{(1/\sigma)+1}}\int_{l^{(2)}}^{k}\frac{(\mu^{(1)}+R^{(1)}\exp(-\sigma(l^{(1)}-k)))^{(1/\sigma)-1}}{(\mu^{(2)}+R^{(2)}\exp(-\sigma(k-l^{(2)})))^{-1/\sigma}}dk$
(37) $\displaystyle+$ $\displaystyle
c^{(2)}(l^{(2)})e^{l^{(2)}-k}\left(\frac{\mu^{(2)}+R^{(2)}}{\mu^{(2)}+R^{(2)}e^{-\sigma(k-l^{(2)})}}\right)^{(1/\sigma)+1}$

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

$\displaystyle c^{(2)}(k)$ $\displaystyle=$ $\displaystyle
c^{(1)}(l^{(1)})\mu^{(1)}e^{-(k-l^{(1)})}\left(\frac{\mu^{(1)}+R^{(1)}}{\mu^{(2)}+R^{(2)}\exp(-\sigma(k-l^{(2)}))}\right)^{(1/\sigma)+1}$
(38)
$\displaystyle\int_{l^{(2)}}^{k}\frac{(\mu^{(2)}+R^{(2)}\exp(-\sigma(k-l^{(2)})))^{1/\sigma}}{(\mu^{(1)}+R^{(1)}\exp(-\sigma(k-l^{(1)})))^{(1/\sigma)+1}}dk$
$\displaystyle+$ $\displaystyle
c^{(2)}(l^{(1)})e^{-(k-l^{(1)})}\left(\frac{\mu^{(2)}+R^{(2)}e^{-\sigma(l^{(1)}-l^{(2)})}}{\mu^{(2)}+R^{(2)}e^{-\sigma(k-l^{(2)})}}\right)^{(1/\sigma)+1}$

The integrals in (34) to (38) can evaluated numerically. The result of the
numerical computation is shown in Fig. A2.

(a)

(b)

Figure A2: Glycan concentration profile calculated from the model using (a)
formula (33) for $N_{E}=N_{C}=1$ and (b) formulae (34)-(38) for
$N_{E}=1,N_{C}=2$.

##  Appendix E Capability of the chemical network model to generate complex
distributions

Is our glycan synthesis model capable of generating concentration
distributions of arbitrary complexity? In what way do we need to change the
parameters $N_{E},N_{C},\sigma,\ldots$, in order to generate glycan
distributions $\mathbf{c}$ of a given complexity? The purpose of this section,
is to obtain some heuristics for this task.

We show in Appendix D that (19) can be solved analytically in the limit
$N_{s}\gg 1$, because in this limit the glycan index $k$ can be approximated
by a continuous variable, and the recursion relations for the steady state
glycan concentrations (5)-(6) can be cast as a matrix differential equation.
This allows us to obtain an _explicit_ expression for the steady state
concentration in terms of the parameters $(\mathbf{R},\mathbf{L})$.

We derive our heuristics from a semi-analytical treatment in the limit
$N_{s}\gg 1$ (Appendix D), which apart from being simple to implement in
general, provides an explicit formula for $c_{k}$ for the case $N_{E}=N_{C}=1$
(33). Figures A3(a)-(d) show the glycan profile $c_{k}\,\,vs.\,k$ as one
varies the enzyme specificity $\sigma$, the reaction rates $R$ and transport
rates $\mu$, for two different values of $N_{E}$ and $N_{C}$. The results in
the plots lead us to the following general observations:

  * •

Very low specificity enzymes cannot generate complex glycan distributions.
Keeping everything else fixed, intermediate or high specificity enzymes can
generate glycan distributions of higher complexity by increasing $N_{E}$ or
$N_{C}$ (Figs. A3(a),(c)).

  * •

Decreasing the specificity $\sigma$ or increasing the rates $R$ increases the
proportion of higher index glycans. Keeping everything else fixed, changes in
the rate $R$ have a stronger impact on the relative weights of the higher
index glycans to lower index glycans. The relative weight of the higher index
glycans increases with increasing $N_{E}$ and $N_{C}$ (Figs. A3(b)-(d)).

  * •

Keeping everything else fixed, decreasing enzyme specificity increases the
spread of the distribution around the peaks (Figs. A3(a),(c)).

(a)

(b)

(c)

(d)

Figure A3: Glycan profile $\\{c_{k}:k=1,\ldots,N_{s}\\}$ as a function of
specificity $\sigma$ (Fig. (a), (c)), and reaction rates $R$ (Fig. (b), (d)).
Fig. (a): $N_{E}=N_{C}=1,(R=50,\mu=1,l=10)$. $c_{k}$ decreases exponentially
with $k$ for very low and very high $\sigma$; however, the decay rate is lower
at low $\sigma$. For intermediate values of $\sigma$, 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
$\sigma$.
Fig. (b): $N_{E}=N_{C}=1,(\sigma=0.1,\mu=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$.
Fig. (c):
$N_{E}=N_{C}=2,(R=40,\mu=1,[l^{(1)}_{1},l^{(1)}_{2},l^{(2)}_{1},l^{(2)}_{2}]=[10,30,50,70])$.
As $\sigma$ increases, the distribution becomes more complex – from a single
peaked distribution at low $\sigma$ to a maximum of four-peaked distribution
at high $\sigma$. The peaks gets sharper, and more well defined as $\sigma$
increases.
Fig. (d):
$N_{E}=N_{C}=2,(R=40,\mu=1,[l^{(1)}_{1},l^{(1)}_{2},l^{(2)}_{1},l^{(2)}_{2}]=[10,30,50,70])$.
As in the plots in Fig. (b), increasing $R$ shifts the peaks towards higher
index glycans and the proportion of higher index glycan increases.

##  Appendix F 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, $\mu=.18$/min. We bound the transport rate
for our optimization between 0.01/min and 1/min.

Next, we estimate the range of values for the chemical reaction rates. The
injection rate $q$ is in the range $100-1500$ pmol/$10^{6}$ cell 24 h Umaña
and Bailey (1997); Krambeck  _et al._ (2009). For our calculation we set
$q=387.30$ pmol/$10^{6}$ cells 24 hr = $0.27$ pmol/$10^{6}$ cells min, the
geometric mean of $100$ and $1500$. We set the range for the enzymatic rate
$R$ to be

$R_{\min}=\min_{\alpha}\Big{\\{}\frac{V^{(\alpha)}_{\max}/\nu}{K^{(\alpha)}_{M}+\frac{1}{\nu}\frac{q}{\mu}}\Big{\\}}\leq
R\leq
R_{\max}=\max_{\alpha}\Big{\\{}\frac{V^{(\alpha)}_{\max}/\nu}{K_{M}^{(\alpha)}}\Big{\\}}.$

where $K^{(\alpha)}_{M}$ and $V^{(\alpha)}_{\max}$ denote the Michaelis
constants and $V_{\max}$ of the $\alpha$-th enzyme. The conversion from 1
pmoles/$10^{6}$ cells to concentration can be obtained by taking cisternal
volume ($\nu$) to be 2.5$\mu m^{3}$ Umaña and Bailey (1997); Krambeck  _et
al._ (2009). This gives

$\text{1 pmoles/$10^{6}$ cells}=\frac{10^{-12}\text{moles}}{10^{6}\times
2.5\times 10^{-18}\times 10^{3}\text{litre}}=400\mu M.$ (39)

In Table 1 we report the parameters for the $8$ enzymes taken from Table 3 in
Umaña and Bailey (1997). From these parameters it follows that

$\displaystyle R_{\min}$ $\displaystyle=$
$\displaystyle\min_{\alpha}\Big{\\{}\frac{V^{(\alpha)}_{\max}/\nu}{K^{(\alpha)}_{M}+\frac{1}{\nu}\frac{q}{\mu}}\Big{\\}}$
$\displaystyle=$
$\displaystyle\frac{V^{(7)}_{\max}/\nu}{K_{M}^{(7)}+\frac{1}{\nu}\frac{q}{\mu}}=\frac{.16\times
400\mu M/\text{min}}{3400\mu M+149.4\mu M}=0.018\text{min}^{-1}$
$\displaystyle R_{\max}$ $\displaystyle=$
$\displaystyle\max_{\alpha}\Big{\\{}\frac{V^{(\alpha)}_{\max}/\nu}{K_{M}^{(\alpha)}}\Big{\\}}$
$\displaystyle=$
$\displaystyle\frac{V^{(1)}_{\max}/\nu}{K^{(1)}_{M}}=\frac{5\times 400\mu
M/\text{min}}{100\mu M}=20\text{min}^{-1}$

$\alpha$ | $K^{(\alpha)}_{M}$ | $V^{(\alpha)}_{\max}$
---|---|---
| ($\mu$mol) | (pmol/$10^{6}$ 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
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^{(\alpha)}_{M}$ and $V^{(\alpha)}_{\max}$ denote the Michaelis constant and
$V_{\max}$ of the $\alpha$-th enzyme.

##  Appendix G Constructing target distributions for glycans of a given cell
type

The target 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 ionised 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.

This 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): the 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 Gaussian Mixture Model (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) had mass ranging between
500 to 5000 Daltons with intensity reported at every 0.0153 Daltons. 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. Fig. A4 shows the raw MS data and the
binned distribution.

Let $\bar{I}_{k}$ represents the relative intensity of the $k$-th bin in the
binned MS graph. We generate a sample population of glycans using the MS data
in the following way:

  1. 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}\bar{I}_{k}$. The plots in Fig. A5(a) show that the results are
not sensitive to the specific value of $I_{\max}$.

  2. 2.

Chemical noise: The sample used for MS analysis also contains small amounts of
molecules that are not glycans. These appear as the 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^{-\bar{I}_{k}}$

which adequately suppresses this chemical noise.

  3. 3.

Bootstrapped glycan data: The count $n_{k}$ at the glycan index $k$ is
distributed according to the following distribution:

$n_{k}=\left\\{\begin{array}[]{lll}0&(1-p_{k})&n=0\\\
n&p_{k}e^{-I_{k}}\frac{(I_{k})^{n}}{n!}&n\geq 1.\end{array}\right.$

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.i.d. samples
distributed according to the distribution above. We in Fig. A5 (b) show that
results are insensitive to $N_{i}$.

Next, we interpret the counts as samples from a “spatial” distribution $f$. We
approximate this distribution as a Gaussian mixture, i.e.
$f(x)=\sum_{\ell=1}^{N}\gamma_{\ell}\,\eta(x\mid\mu_{\ell},\sigma_{\ell})$,
where $\eta(x\mid\mu,\sigma)$ denotes the density of a normally distributed
random variable with mean $\mu$ and standard deviation $\sigma$ at the
location $x$. In this setting, we assume that each count is a sample from the
distribution $\eta(x\mid\mu_{\ell},\sigma_{\ell})$ with probability
$\gamma_{\ell}$. Thus, each count is classified as coming from one of the
Gaussian components.

Figure A4: The binned MS data (blue) approximates the raw MS data (red) very
well. We use this binned data for GMM approximation of the MS data.

(a)

(b)

Figure A5: Log likelihood vs. number of components ($N$) in the GMM. We see
that the log likelihood saturates at around $N=20$, thus $20$-GMM is a very
good representation of the MS-data from human T-cells. The different symbols
are for (a) different values of the maximum intensity $I_{max}=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_{max}$ and
$N_{i}$.

##  Appendix H Numerical scheme for performing the non-convex optimization

We solve Optimization C using the numerical scheme detailed below. The
optimization problem consists of minimising a non-convex objective with linear
box constraints. We use the MATLAB FMINCON function to solve this
optimization. We use Sequential Quadratic Programming (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 initialise the algorithm with many random
initial points. We use SOBOLSET function of MATLAB to generate space filling
pseudo random numbers. We have taken $1000$ initialisations for each
$N_{E},N_{C}$ and $\sigma$ value. We have taken $50$ equally spaced points
between $0$ and $1$ to explore the $\sigma$-space for Fig. 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.
$\sigma$ graph. The results for $\sigma_{min}(N_{E},N_{C})$ and
$D(\sigma_{min},N_{E},N_{C})$ (Fig. 4) were obtained by adding $\sigma$ to the
optimization vector and then performing the optimization. The sensitivity
results (Figs. 4e and 4e) were obtained by approximating the $D$ vs $\sigma$
graph around $\sigma_{min}$ with a parabola, the coefficient of the quadratic
term being the curvature of the graph at $\sigma_{min}$.

A similar numerical scheme was used to optimize diversity.

(a) $c_{th}=\frac{1}{N_{s}}$

(b) $c_{th}=\frac{1}{2N_{s}}$

(c) $c_{th}=\frac{1}{4N_{s}}$

Figure A6: Diversity vs. $N_{C}$ for different values of $\sigma$ keeping
$N_{E}=1$ fixed, for three different values of the threshold,
$c_{th}=\frac{1}{N_{s}},\;\frac{1}{2N_{s}},\;\frac{1}{4N_{s}}$. Changing the
value of the threshold $c_{th}$, only changes the saturation value of the
diversity curve.

## References

  * Alberts  _et al._ (2002) B. Alberts _et al._ ,  _Molecular Biology of the Cell_ (Garland Science, 2002).
  * Varki  _et al._ (2009) A. Varki _et al._ ,  _Essentials of Glycobiology_ (Cold Spring Harbor Laboratory Press, 2009).
  * Cummings and Pierce (2014) R. D. Cummings and J. M. Pierce, Chemistry & biology 21, 1 (2014).
  * Varki (2017) A. Varki, Glycobiology 27, 3 (2017).
  * Drickamer and Taylor (1998) K. Drickamer and M. E. Taylor, Trends in biochemical sciences 23, 321 (1998).
  * Gagneux and Varki (1999) P. Gagneux and A. Varki, Glycobiology 9, 747 (1999).
  * Gabius (2018) H.-J. Gabius, BioSystems 164, 102 (2018).
  * Dwek (1996) R. A. Dwek, Chemical reviews 96, 683 (1996).
  * Winterburn and Phelps (1972) P. Winterburn and C. Phelps, Nature 236, 147 (1972).
  * Varki (1998) A. Varki, Trends in cell biology 8, 34 (1998).
  * Pothukuchi  _et al._ (2019) P. Pothukuchi, I. Agliarulo, D. Russo, R. Rizzo, F. Russo, and S. Parashuraman, FEBS letters 593, 2390 (2019).
  * Bard and Chia (2016) F. Bard and J. Chia, Trends in cell biology 26, 379 (2016).
  * D’Angelo  _et al._ (2013) G. D’Angelo, S. Capasso, L. Sticco, and D. Russo, The FEBS journal 280, 6338 (2013).
  * Demetriou  _et al._ (2001) M. Demetriou, M. Granovsky, S. Quaggin, and J. W. Dennis, Nature 409, 733 (2001).
  * WILLS and GREEN (1995) C. WILLS and D. R. GREEN, Immunological reviews 143, 263 (1995).
  * Cover and Thomas (2012) T. M. Cover and J. A. Thomas,  _Elements of information theory_ (John Wiley & Sons, 2012).
  * MacKay (2003) D. J. MacKay,  _Information theory, inference and learning algorithms_ (Cambridge university press, 2003).
  * Sengupta and Linstedt (2011) D. Sengupta and A. D. Linstedt, Annual review of cell and developmental biology 27, 57 (2011).
  * Sachdeva  _et al._ (2016) H. Sachdeva, M. Barma, and M. Rao, Scientific reports 6, 1 (2016).
  * Sens and Rao (2013) P. Sens and M. Rao, in  _Methods in cell biology_ , Vol. 118 (Elsevier, 2013) pp. 299–310.
  * Bacharoglou (2010) A. Bacharoglou, Proceedings of the American Mathematical Society 138, 2619 (2010).
  * Cummings and Crocker (2020) R. D. Cummings and P. Crocker,  _Functional Glycomics Database, Consortium for Functional Glycomics,_ http://www.functionalglycomics.org (2020).
  * Umaña and Bailey (1997) P. Umaña and J. E. Bailey, Biotechnology and bioengineering 55, 890 (1997).
  * Krambeck  _et al._ (2009) F. J. Krambeck, S. V. Bennun, S. Narang, S. Choi, K. J. Yarema, and M. J. Betenbaugh, Glycobiology 19, 1163 (2009).
  * Krambeck and Betenbaugh (2005) F. J. Krambeck and M. J. Betenbaugh, Biotechnology and Bioengineering 92, 711 (2005).
  * Fisher  _et al._ (2019) P. Fisher, H. Spencer, J. Thomas-Oates, A. J. Wood, and D. Ungar, Cell reports 27, 1231 (2019).
  * Fisher and Ungar (2016) P. Fisher and D. Ungar, Frontiers in cell and developmental biology 4, 15 (2016).
  * Hirschberg  _et al._ (1998) C. B. Hirschberg, P. W. Robbins, and C. Abeijon,  _Transporters of nucleotide sugars, ATP, and nucleotide sulfate in the endoplasmic reticulum and Golgi apparatus,_ (1998).
  * Caffaro and Hirschberg (2006) C. E. Caffaro and C. B. Hirschberg, Accounts of chemical research 39, 805 (2006).
  * Berninsone and Hirschberg (2000) P. M. Berninsone and C. B. Hirschberg, Current opinion in structural biology 10, 542 (2000).
  * Trinajstic (2018) N. Trinajstic,  _Chemical graph theory_ (Routledge, 2018).
  * Price and Stevens (1999) N. Price and L. Stevens,  _Fundamentals of Enzymology: The cell and molecular biology of catalytic proteins_ (Oxford University Press, 1999).
  * Moremen and Haltiwanger (2019) K. W. Moremen and R. S. Haltiwanger, Nature chemical biology 15, 853 (2019).
  * Kellokumpu (2019) S. Kellokumpu, Frontiers in cell and developmental biology 7, 93 (2019).
  * Casey  _et al._ (2010) J. R. Casey, S. Grinstein, and J. Orlowski, Nature reviews Molecular cell biology 11, 50 (2010).
  * Dmitrieff  _et al._ (2013) S. Dmitrieff, M. Rao, and P. Sens, Proceedings of the National Academy of Sciences 110, 15692 (2013).
  * Llopis  _et al._ (1998) J. Llopis, J. M. McCaffery, A. Miyawaki, M. G. Farquhar, and R. Y. Tsien, Proceedings of the National Academy of Sciences 95, 6803 (1998).
  * Monod  _et al._ (1965) J. Monod, J. Wyman, and J.-P. Changeux, J Mol Biol 12, 88 (1965).
  * Changeux and Edelstein (2005) J.-P. Changeux and S. J. Edelstein, Science 308, 1424 (2005).
  * Savir and Tlusty (2007) Y. Savir and T. Tlusty, PloS one 2, e468 (2007).
  * Roseman (2001) S. Roseman, Journal of Biological Chemistry 276, 41527 (2001).
  * Hossler  _et al._ (2007) P. Hossler, B. C. Mulukutla, and W.-S. Hu, PloS one 2 (2007).
  * Yang  _et al._ (2018) M. Yang, C. Fehl, K. V. Lees, E.-K. Lim, W. A. Offen, G. J. Davies, D. J. Bowles, M. G. Davidson, S. J. Roberts, and B. G. Davis, Nature chemical biology 14, 1109 (2018).
  * Bar-Even  _et al._ (2015) A. Bar-Even, R. Milo, E. Noor, and D. S. Tawfik, Biochemistry 54, 4969 (2015).
  * Boyd and Vandenberghe (2004) S. Boyd and L. Vandenberghe,  _Convex optimization_ (Cambridge university press, 2004).
  * Jaiman and Thattai (2018) A. Jaiman and M. Thattai, BioRxiv , 440792 (2018).
  * Du  _et al._ (2008) P. Du, G. Stolovitzky, P. Horvatovich, R. Bischoff, J. Lim, and F. Suits, Bioinformatics 24, 1070 (2008).
  * Varki (2011) A. Varki, Cold Spring Harbor perspectives in biology 3, a005462 (2011).
  * Dennis  _et al._ (2009) J. W. Dennis, I. R. Nabi, and M. Demetriou, Cell 139, 1229 (2009).
  * van Halbeek  _et al._ (1983) H. van Halbeek, G. J. Gerwig, J. F. Vliegenthart, H. L. Smits, P. J. Van Kerkhof, and M. F. Kramer, Biochimica et Biophysica Acta (BBA)-Protein Structure and Molecular Enzymology 747, 107 (1983).
  * McFarlane  _et al._ (2014) H. E. McFarlane, A. Döring, and S. Persson, Annual review of plant biology 65, 69 (2014).
  * Koch  _et al._ (2015) B. E. Koch, J. Stougaard, and H. P. Spaink, Glycobiology 25, 469 (2015).
  * O’Neill  _et al._ (2004) M. A. O’Neill, T. Ishii, P. Albersheim, and A. G. Darvill, Annu. Rev. Plant Biol. 55, 109 (2004).
  * Hayashi and Kaida (2011) T. Hayashi and R. Kaida, Molecular Plant 4, 17 (2011).
  * Kumar  _et al._ (2011) P. Kumar, M. Yang, B. C. Haynes, M. L. Skowyra, and T. L. Doering, Current opinion in structural biology 21, 597 (2011).
  * Gow and Hube (2012) N. A. Gow and B. Hube, Current opinion in microbiology 15, 406 (2012).
  * Atmodjo  _et al._ (2013) M. A. Atmodjo, Z. Hao, and D. Mohnen, Annual review of plant biology 64 (2013).
  * Free (2013) S. J. Free, in  _Advances in genetics_ , Vol. 81 (Elsevier, 2013) pp. 33–82.
  * Pauly  _et al._ (2013) M. Pauly, S. Gille, L. Liu, N. Mansoori, A. de Souza, A. Schultink, and G. Xiong, Planta 238, 627 (2013).
  * Burton and Fincher (2014) R. A. Burton and G. B. Fincher, Frontiers in plant science 5, 456 (2014).
  * Becker and Melkonian (1996) B. Becker and M. Melkonian, Microbiol. Mol. Biol. Rev. 60, 697 (1996).
  * Mironov  _et al._ (2017) A. A. Mironov, I. S. Sesorova, E. V. Seliverstova, and G. V. Beznoussenko, Tissue and Cell 49, 186 (2017).
  * Donohoe  _et al._ (2007) B. S. Donohoe, B.-H. Kang, and L. A. Staehelin, Proceedings of the National Academy of Sciences 104, 163 (2007).
  * Mogelsvang  _et al._ (2003) S. Mogelsvang, N. Gomez-Ospina, J. Soderholm, B. S. Glick, and L. A. Staehelin, Molecular biology of the cell 14, 2277 (2003).
  * Ladinsky  _et al._ (2002) M. S. Ladinsky, C. C. Wu, S. McIntosh, J. R. McIntosh, and K. E. Howell, Molecular biology of the cell 13, 2810 (2002).
  * Stanley (2011) P. Stanley, Cold Spring Harbor perspectives in biology 3, a005199 (2011).
  * Nam  _et al._ (2012) H. Nam, N. E. Lewis, J. A. Lerman, D.-H. Lee, R. L. Chang, D. Kim, and B. O. Palsson, Science 337, 1101 (2012).
  * Peracchi (2018) A. Peracchi, Trends in biochemical sciences (2018).
  * Blanes  _et al._ (2009) S. Blanes, F. Casas, J. Oteo, and J. Ros, Physics Reports 470, 151 (2009).