Stochastic modelling of a three-dimensional glycogen granule synthesis and impact of the branching enzyme

In humans, glycogen storage diseases result from metabolic inborn errors, and can lead to severe phenotypes and lethal conditions. Besides these rare diseases, glycogen is also associated to widely spread societal burdens such as diabetes. Glycogen is a branched glucose polymer synthesised and degraded by a complex set of enzymes. Over the past 50 years, the structure of glycogen has been intensively investigated. Yet, the interplay between glycogen structure and the related enzymes is still to be characterised. In this article, we develop a stochastic coarse-grained and spatially resolved model of branched polymer biosynthesis following a Gillespie algorithm. Our study largely focusses on the role of the branching enzyme, and first investigates the properties of the model with generic parameters, before comparing it to in vivo experimental data in mice. It arises that the ratio of glycogen synthase over branching enzyme activities drastically impacts the structure of the granule. We deeply investigate the mechanism of branching and parametrise it using distinct lengths. Not only do we consider various possible sets of values for these lengths, but also distinct rules to apply them. We show how combining them finely tunes glycogen macromolecular structure. Comparing the model with experimental data confirms that we can accurately reproduce glycogen chain length distributions in wild type mice. Additional granule properties obtained for this fit are also in good agreement with typically reported values in the experimental literature. Nonetheless, we find that the mechanism of branching must be more flexible than usually reported. Overall, we demonstrate that the chain length distribution is an imprint of the branching activity and mechanism. Our generic model and methods can be applied to any glycogen data set, and could in particular contribute to characterise the mechanisms responsible for glycogen storage disorders. Author summary Glycogen is a granule-like macromolecule made of 10,000 to 50,000 glucose units arranged in linear and branched chains. It serves as energy storage in many species, including humans. Depending on physiological conditions (hormone concentrations, glucose level, etc.) glycogen granules are either synthesised or degraded. Certain metabolic disorders are associated to abnormal glycogen structures, and structural properties of glycogen might impact the dynamics of glucose release and storage. To capture the complex interplay between this dynamics and glycogen structural properties, we propose a computational model relying on the random nature of biochemical reactions. The granule is represented in three dimensions and resolved at the glucose scale. Granules are produced under the action of a complex set of enzymes, and we mostly focus on those responsible for the formation of new branches. Specifically, we study the impact of their molecular action on the granule structure. With this model, we are able to reproduce structural properties observed under certain in-vivo conditions. Our biophysical and computational approach complements experimental studies and may contribute to characterise processes responsible for glycogen related disorders.


Introduction
1 Management of energy resources is crucial for an organism to survive, since nutrient 2 availability may vary considerably with time. Moreover, organisms face numerous 3 other stresses that may temporarily increase energy demand. To react to such changes 4 in energy supply and demand, it is apparent that some internal stores are necessary. 5 These stores are filled when nutrients are abundant and depleted when demand 6 exceeds available supply. As direct substrate for catabolic pathways, glucose plays a 7 central role in energy metabolism in most organisms [1]. While plants have found their 8 way to store glucose as starch, animals, fungi, and most bacteria store glucose as 9 glycogen. Both starch and glycogen are branched polymers consisting of glucose 10 monomers, linked into linear chains by α-1,4 bonds, and branching points by α-1,6 11 bonds. However, these two macromolecules exhibit rather different physico-chemical 12 properties. In contrast to glycogen, starch is insoluble in water under physiological 13 conditions, and contains high density crystalline regions. These different properties are 14 reflected by distinct branching patterns and chain length distributions. Functionally, 15 starch and glycogen can be compared to capacitors in electric circuits. The latter are 16 able to store and release electrons depending on current and voltage. Thus, they can 17 be used to stabilise a fluctuating electric signal. Analogously, glycogen and starch can 18 be seen as two different capacitors that both contribute to glucose homeostasis by 19 managing energy resources through time. 20 While the fine structure of starch has been widely investigated over the past 50 21 years, less is known on that of glycogen [2]. The characterisation of the detailed 22 structure of glycogen, as well as the interplay between its structural properties and the 23 dynamics of glycogenesis (synthesis) and glycogenolysis (degradation) is unclear. Yet, 24 both glycogen structure and dynamics are of utmost interest for understanding 25 glycogen metabolism and the impact of related genetic variances. For human health, 26 this is particularly important considering the increasing prevalence of glycogen storage 27 diseases (GSDs), as well as diabetes, and other glycogen related disorders. whose synthesis and degradation both involve a relatively small number of enzymes 41 (see Fig 1). For a comprehensive review, see [16]. 42 To initiate a new glycogen molecule, a chain of 5 glucose units is synthesised and 43 bound to a glycogenin (GN) protein, which forms the core of the final granule [17][18][19][20]. 44 Once initiation is completed, the granule is expanded by the two enzymes Glycogen 45 Synthase (GS) and Glycogen Branching Enzyme (GBE). GS is an elongating enzyme 46 that adds one glucose residue to the non-reducing end of an α-1,4 linear chain, thereby 47 forming a new α-1,4 glucosydic bond. GBE cleaves a small part of a linear chain and 48 creates a new branch by forming an α-1,6 glucosydic bond. We call "daughter" the 49 newly formed chain and "mother" the one it is branching from. 50 Besides synthesis, granules are subject to degradation, that is performed by two 51 other enzymes. Glycogen Phosphorylase (GP) and Glycogen Debranching Enzyme 52 (GDE) respectively shorten and debranch glycogen branches. 53 Depending on the relative enzymatic activities of these four enzymes (GS, GBE, 54 GP, and GDE) the overall size of the glycogen granule can either increase or decrease. 55 In this article, we choose to focus on glycogen synthesis, and more specifically the role 56 of the branching enzyme GBE. In vivo, the GS and GBE enzymes synthesise glycogen, while the GP and GDE degrade it. Besides, GN is the initial precursor of the granule and stands in its core. Enzymes are noted in orange, glucose residues are in blue, and GN is highlighted with a yellow sphere.
Experimental observations [21] of glycogen show an average chain length (CL) of 58 13 glucose units. The chain length distribution is not symmetric [22,23]. A typical 59 peak is observed at low degree of polymerisation (DP), around DP 8, and almost no 60 chains are detected above DP 40. The degree of branching is defined in two ways in 61 the literature. It is most commonly defined as the ratio of α-1,6 to α-1,4 linkages, but 62 sometimes also as the average number of α-1,6 bonds per chain [6]. We will apply the 63 first definition throughout the paper. This ratio is in the range 0.02 − 0.05 in 64 amylopectin [24][25][26] from starch and 0.06 − 0.10 in glycogen [9,27]. In 1956, Peat et explaining macroscopic features of glycogen granules. One of them has become known 72 as the Whelan model, which assumes that every chain carries on average two branches. 73 Based on the Whelan model, the nowadays widespread idea emerged that glycogen 74 can be described as a fractal molecule [6,31,32]. A fractal glycogen structure is indeed 75 attractive because it can reproduce various structural properties of glycogen. 76 Moreover, it provides a mechanism explaining why glycogen granules seem to have a 77 maximum size. With this model, glycogen becomes too dense at the surface due to an 78 exponential increase of the number of non-reducing ends with the distance from the 79 centre. Thus, steric hindrance prevents synthesis to continue. However, Manners [2] 80 stressed in 1990 that there is no clear evidence that supports a regular branching and 81 therefore a fractal pattern. More recently, further arguments and results against a 82 fractal structure have been raised [33,34]. Besides, independently from any fractal 83 structure considerations, Zhang et al. [35]  Glycogen structure and geometry 101 We represent glycogen granules as simplified three-dimensional structures, in which 102 every glucose monomer is characterised by its position in space (see Fig 2). To describe 103 the branched tree-like structure, we generalise the simple representation of linear 104 self-avoiding polymers. Using X-ray experiments, Goldsmith et al [31] characterised in 105 detail how glucose molecules are arranged into helical α-1,4 linear chains. The 106 cross-section of the helix has been calculated to be 1.3 nm 2 with 6 to 7 residues per 107 turn of length 1.4 nm. The radius of the circular cross-section is thus ρ = 0.65 nm, 108 and each glucose residue contributes to the chain length by l = 0.24 nm. Inspired by 109 these findings, we propose that monomers are described as overlapping spheres with 110 radius ρ = 0.65 nm, equal to that of the helix. The validity of this assumption and its 111 impact on our results are presented in detail in the Supplementary Material (see Fig   112   11). Besides, the center of consecutive monomers are distant by l = 0.24 nm, which 113 corresponds to the contribution of one glucose unit to the chain length, but also 114 involves that the coarse-grained monomer spheres overlap. However, to realistically 115 reflect self-avoidance, two different chains cannot. With these spatial considerations, 116 we ensure that the contribution to the volume by one glucose unit in the model is 117 similar to that of the real helical chain. This way, we provide a description which is 118 simple enough to be easily represented in a computer model, but still realistic enough 119 to reflect the spatial properties of linear chains arranged into helices. Assuming helical chains, glucose units are described as interpenetrated spheres with radius ρ = 0.65 nm. Two consecutive glucoses are distant by l = 0.24 nm, which is the radial contribution to the chain length of one glucose in a helical structure. Right: Description of a branching point. We generate the direction of the new branch by randomly picking two angles φ and ψ. The first monomer of the new branch will be located at a distance greater than 2ρ to insure no overlapping between the mother and the daughter branches.
Describing branches formed by α-1,6 linkages requires additional geometrical 121 considerations. As illustrated in Fig 2 (right panel), a branch point is defined by the 122 monomer on the mother chain to which it is attached, and two angles defining the 123 direction of the daughter chain. Besides, the anchoring monomer on the mother chain 124 and the first one of the daughter chain are distant by 2ρ, ensuring that they do not 125 overlap.

126
Glycogen synthesis is initiated by GN that is located in the centre of the 127 molecule [36]. It contributes to the total volume of glycogen, and thus we also consider 128 self-exclusion between GN and any glucose residue of the granule. For the sake of 129 simplicity, we assume that it is a sphere of density 1.3 g · cm −3 [37], which is the 130 typical density of a protein. Accounting for its two sub-units of 38 kDa each [11], we 131 approximate its radius by ρ GN ≈ 2.85 nm. Together with GS, which is responsible for 132 elongation, GN may initiate more than a single primary chain, possibly either 2 or 4 133 [20]. We model this initial core structure with two primary chains pointing out of the 134 GN sphere in opposite directions.  GS binds the non-reducing end of an α-1,4 linear chain and elongates it by adding 139 one glucose residue. It is commonly assumed that GS needs a glucose chain primer as 140 a glucose acceptor [38]. In the model, we call L GS min the minimal chain length of this 141 required primer by GS. As elongation takes place, the chain becomes long enough to 142 be cleaved and branched. This reaction is catalysed by GBE. Just like GS, the action 143 of GBE is also characterised by specific substrate and product lengths as shown in Fig 144  3. Since less is known for GBE, we tested two models for its mechanism, namely the process, and creating a new A chain. To precisely describe this last step, we define two 156 additional lengths. First, the new α-1,6 branching point must not be closer than 157 L GBE spacing from either the first glucose of the chain, or an above α-1,6 branching point.

158
Second, the new branching point must not be closer than L GBE leftover from the 159 non-reducing end of the substrate chain, which is the original position of cleavage.

160
Thus, for a new branch to be created by GBE, the substrate branch must have a chain 161 length greater than L GBE min , verifying: To illustrate the impact of these minimal lengths on the reaction outcomes, in Fig   163   4 we detail the case of {L GBE spacing , L GBE transferred , L GBE leftover } = {2, 2, 2}. If the substrate 164 reaches a length of DP equal 7, the condition (Eq 1) is fulfilled and the reaction may 165 take place. If this reaction occurs, there is a single possible outcome (Fig 4, left panel). 166 If the branching reaction occurs on a longer chain than just the minimal one, several  With these minimal lengths, the minimal DP required for a branching to occur is DP = 7. If the chain length is longer, the number of possible outcomes increases. Left: With a substrate of DP = 7, only one outcome is possible. Right: With a substrate of DP = 9, up to 6 distinct outcomes are possible.

172
When the simulation starts, the system is composed of a GN core with two primary 173 chains standing in opposite directions in the center of the simulation space. Two 174 enzymes (GS and GBE) modify the structure of the glycogen granule. To quantify 175 that, we define Γ as the ratio of the GS over GBE reaction rates. For Γ ≈ 0, branching 176 dominates over elongation, and vice versa when Γ ≫ 1.

Granule density 197
Our approach tracks the position (x, y, z) of each glucose unit in three-dimensions. This allows us to compute how densely the granules are packed. Granule packing is quantified by the occupancy Ω, which is defined by the volume occupied by glucoses V glucose divided by the total volume V total , To characterise an entire granule, we consider the total volume V total to be a sphere of 198 gyration radius R g (defined in the section 5.1 of the Supplementary Material). To 199 determine the occupancy at a given radius r from the center of the granule, we 200 estimate the local occupancy in a spherical shell between radii r and r + ∆r.

216
The relatively low occupancy close to the granule center is due to the presence of 217 the GN core, which is not counted in the occupancy, but its corresponding volume 218 cannot be filled with glucose units. We observe that at most 30% of the granule 219 volume is occupied by glucose. This value is rather low, which is expected, since in the 220 model, branches are straight helices, without any flexibility, while in reality, branches 221 can bend and form locally higher densities. Specifically, this occupancy value is below 222 Ω = 0.45, which we can estimate from previous studies [31]. There, the authors 223 combine experimental and modelling approaches to conclude that granules of 55,000 224 glucoses have a radius of 21 nm [31]. Nonetheless, there are many uncertainties on the 225 molecular masses experimentally measured. Thus, the occupancy value (i.e. Ω = 0.45) 226 that we deduce from their work might hold large errors. For instance, it is unclear how 227 the water and protein molecules embedded in the granules contribute to the molecular 228 masses experimentally measured.

230
The branching enzyme mechanism is characterised by a triplet of integer numbers, 231 denoted {L GBE spacing , L GBE transferred , L GBE leftover }, which specifies a unique set of rules for the 232 enzymatic reaction. These rules considerably impact the CLD. Fig 7 shows that when 233 these minimal lengths increase, the peak of the CLD is less pronounced and the 234 distribution spreads towards higher DPs. Also, a change in each minimal length has a 235 specific effect on the CLD.

236
Increasing L GBE spacing drastically reduces the peak and spreads the distribution, while  It is important to notice, that these results are obtained when branching dominates 246 over elongation. Increasing the elongation to branching ratio Γ systematically 247 smoothens any multi-modal CLD, because it introduces flexibility in the branching 248 location. It also flattens the peak and spreads the distribution towards higher DPs.

249
Consequently, bi-modal distributions are obtained for high L GBE spacing and low Γ values.

250
As discussed in the Elongation to branching ratio section, we can compare the 251 synthesis process of glycogen and starch. Specifically, in starch the CLD is bi-modal 252 and shifted towards higher DPs as compared to glycogen [40]. Based on our preceding 253 remarks, it could mean that starch branching enzymes are characterised by large 254 substrate specificity lengths, corresponding to a more constrained mechanism than for 255 glycogen. An alternative explanation for the arising of multi-modal CLDs, also based  , respectively. When L GBE spacing increases, a multi-modal distribution emerges. Top-right: CLD for L GBE leftover = 2, 4, and 6, respectively. When L GBE leftover increases, the peak is reduced and the overall distribution spreads towards higher DPs. Bottom-left: CLD for L GBE transferred = 2, 4, and 6, respectively. When L GBE transferred increases, the distribution shifts towards higher DPs. Bottom-right: CLD for L GBE spacing = L GBE transferred = L GBE leftover = 2, 4, and 6, respectively. Varying these distinct minimal lengths concomitantly, combines the individual effects described above, when a single length is varied. Each CLD is the result of averaging 200 simulations of granules with 5,000 glucose units each.
that we compare to experimental results.

268
GBE's mechanism and substrate specificity are incompletely characterised, yet they 269 drastically impact on the CLD, as shown in Fig 7. Therefore, we specifically focus on 270 this enzyme. To this end, we use experimental data obtained by Sullivan et al. [23] for 271 mice muscles, that we extracted from their publication using the software Engauge 272 Digitizer [41]. After purification, the granule chains are debranched using isoamylase, 273 and their degree of polymerisation is measured using size exclusion chromatography.

274
Our fitting procedure can be applied to any glycogen data, obtained from any specie 275 and tissue. The data by Sullivan et al. [23] present two major advantages for our study. 276 First, they are quantitative measurements of good resolution. Second, the authors 277 investigated the case of a defective GBE.
spacing is represented on the Y-axis, ranging from 1 to 6. The elongation to branching ratio Γ is varied on the X-axis. We notice that only low values for L GBE spacing , L GBE leftover and L GBE transferred allow to fit the experimental data. However, one should remember that these are minimal lengths and that the positions at which GBE is able to cleave and branch is flexible beyond these minimal lengths (see Fig 4).
Overall, not only have we ran our fitting procedure using both the flexible and the strict location branching models, but also considering distinct values for ρ (see Figs 14 and 15 in the Supplementary Material). Interestingly, the best fit is obtained for ρ = 0.65 nm, which reflects the realistic size of a glucose unit inside an helical chain. This highlights the role of the steric hindrance to mimic the granule's structural properties. The best fit is obtained for the following set of parameters: Importantly, the parameter values for the branching enzyme inferred from the best fit 279 are distinct from those reported in the field [31,42,43], especially for the typical 280 spacing observed between two branches [31]. Also, based on our results, GBE is able 281 to transfer less than 4 glucose units. Knowing that GS's chain length specificity has to 282 follow the same rule, this questions the commonly assumed value of DP4 as the 283 minimal length that can be elongated by GS.

284
In Fig 9,  The Y -axis shows L GBE spacing ranging from 1 to 6. The X-axis shows the elongation to branching ratio Γ ranging from 0.1 to 5.0. A given cell corresponds to a parameter set {Γ, L GBE spacing , L GBE leftover , L GBE transferred , ρ}. Additional parameter sets are tested around good scores, i.e. the resolution on the elongation to branching ratio is increased, as well as the number of runs averaged. This area is surrounded by a red rectangle in which the average score is 12.26. Fitting scores are ranging from 8.6 to 369.0. The best score is 8.6 (red square in the inset heat-map) which corresponds to {Γ = 0.6, L GBE spacing = 1, L GBE transferred = L GBE leftover = 3, ρ = 0.65 nm}.

313
It is interesting to notice that ca. 21 nm is the upper limit for the radius of a 314 glycogen granule of ca. 50,000 glucoses, when considering the fractal structure 315 depicted by Meléndez et al. [6]. Thus, with our approach, that is not based on a 316 fractal structure, we determine a gyration radius close to their results. However, the 317 glycogen structures we simulate are deeply distinct from theirs. In our simulations, we 318 observe that the density is approximately constant with the radius, while a fractal 319 glycogen granule has a density that increases exponentially with the radius, so 320 resembling a rather empty shell. Also, we show that our simulations can reproduce 321 characteristic quantities (summarised in Table 1 for Γ = 0.6), which are in good 322 agreement with experimental results. Thus, our model appears to provide a more 323 realistic depiction of glycogen granules.

324
For comparison, in Table 1, we also report the case of Γ = 50.0, while keeping all 325 other parameters unchanged. By increasing Γ, we force elongation over branching. As 326 expected, less chains are created but they are longer. Consistently, the branching degree is also lower. In our model, chains are rigid rods that cannot bend, which leads 328 to a bigger radius of gyration. Since the total content of glucose is the same like for 329 Γ = 0.6, the occupancy Ω decreases almost proportionally. In vivo, glycogen granules 330 might follow a distinct behaviour, where the presence of long chains could trigger 331 precipitation, like for instance reported in the Lafora and GBE diseases [23,45].

332
The formation of an A chain results from branching either an A or a B chain. These different events can be represented as reactions and described as follows: To illustrate how different branch structures correspond to specific A:B ratios, we 333 sketch extreme cases in Fig 10. If the R B→AB and R A→AB reactions are equally likely, 334 a purely probabilistic approach predicts that, on average, the A:B ratio is equal to 1, has a typical A:B ratio in the range from 1.5 to 2.6 [29,30], depending on the 341 organism. It is therefore interesting to investigate which factors influence this ratio.

342
Here, we identify two effects of the interplay between the branching mechanism and 343 the dynamics of the granule structure, highlighting once again the importance of the 344 branching enzyme on the emerging structural patterns. First, considering a small Γ, 345 branching dominates over elongation, and the enzyme branches as soon as it can. 346 Thus, if in addition L GBE spacing is small, the A chains are closely stacked together and the 347 granule is tightly packed. In this regime, overall, A chains are closer to other A chains 348 than to B chains. Due to the increased steric hindrance for A chains, B chains react 349 more easily, which favours the R B→AB reaction, and leaves A chains unchanged.

350
Therefore, the A:B ratio increases, meaning that it shifts to the right along the axis in 351 Fig 10. Second, if L GBE transferred > L GBE leftover , when a branching occurs, on either an A or a 352 B chain, the newly formed A chain will on average be longer than the B (leftover) one. 353 Therefore, the new A chain will react faster, favouring R A→AB , and lead to a  average a single glycogen synthase enzyme is active per granule [46]. As a result, a 378 stochastic approach appears very natural to deal with this complex and spatially 379 structured substrate. 380 We base our stochastic approach on the Gillespie algorithm [47]. At each time step, 381 it consists in randomly selecting both an enzymatic reaction and its duration, while

404
The best fit algorithm consists in setting bounding ranges for the parameters to be fitted, and systematically scanning the parameter space thereby defined. For each parameter set tested, we run 50 simulations, take the average of the resulting CLDs, and compare it to our targeted experimental data set. For each parameter set, the comparison to the experimental data is measured by the score S, defined as follows: where A exp (DP) is the experimental CLD abundance for a given DP andĀ sim (DP) 405 the average simulated CLD abundance for a given DP. Therefore, S measures the 406 difference between the average simulated and the experimental CLDs. The best fit is 407 found for the parameter set that has the lowest score S. it is interesting to note that a glycogen fractal description would instead give a density 441 exponentially increasing with the granule radius. Besides our own results and the 442 various arguments exhibited against a fractal representation of glycogen, we would like 443 to highlight here that, as soon as one considers even a single degree of freedom in the 444 glycogen branching reaction, it would rapidly lead to the loss of any "fractal-like" 445 structure. Such degrees of freedom are necessarily present in natural conditions. For 446 instance, the dihedral angles defining the α − 1, 6 bonds may take various values, 447 making a fractal pattern very unlikely in vivo. Beyond these spatial considerations, the 448 model establishes a natural and clear connection between enzymes' reaction rates, and 449 both the degree of polymerisation of the chains, and the number of non-reducing ends. 450 Our results show that the chain length distribution (CLD) of glycogen is highly 451 sensitive to the branching reaction, predominantly its mechanism. Each of the three 452 characteristic minimal lengths of the reaction has specific effects on the CLD.

453
Additionally, if any of them increases, less branching outcomes are possible, eventually 454 leading to a bi-modal or even multi-modal distribution. This effect is enhanced by 455 high branching activity. When varied together, these minimal lengths show even more 456 complex imprints on the CLD. In contrast, multi-modal distributions become less 457 pronounced by increasing the elongation reaction rate, because rapid elongation 458 increases the number of possible configuration outcomes. Additionally, increasing 459 elongation leads to longer chains and results in a CLD spreading towards higher DPs.

460
Altogether, we show that the CLD, and in particular peaks location and intensity, are 461 subtly affected by several complex effects.

462
Guided by these findings, we propose to consider the CLD not only as an 463 important structural feature of glycogen, but also as a signature of GBE. Thus, fitting 464 experimental data with our model arises as a natural strategy to infer knowledge on 465 the GBE mechanism. Not only did we illustrate the strength of our fitting procedure 466 on an experimental data set from Sullivan et al. [23], but we also extracted several 467 macroscopic characteristic features from the resulting fit, which we compared to 468 various literature sources, thereby confirming the validity of our results. Using this fit, 469 at the microscopic level, we were able to discriminate between the two branching  Moreover, our customisable branching model shows that the A:B ratio is 481 independent of the kinetic parameter Γ, but specifically determined by the difference 482 between L GBE transferred and L GBE leftover . Based on these observations, we hypothesise that the 483 branching mechanism is chiefly responsible for the structural differences observed 484 between starch and glycogen.

485
Throughout this study, our coarse-grained approach accounts for the contribution 486 of individual glucoses to the overall granule structure, by considering them as spheres 487 of radius 0.65 nm. It is interesting to notice that if we set the glucose volume to 0 nm 3 , 488 the CLD remains almost unchanged while other macroscopic properties of the granule, 489 such as its overall volume, are dramatically impacted. This observation confirms once 490 more that the glycogen CLD is primarily shaped by the enzyme mechanisms.

492
It is important to keep in mind that our model contains limitations that may be 493 circumvented by further refining the model assumptions. In our model, the simplified 494 coarse-grained representation of glucoses assumes that all of them are arranged in 495 single helices. This hypothesis implies that all α − 1, 4 glycosidic bonds have the same 496 angle values. Yet, in vivo, this is highly unlikely, and instead, the dihedral angles of 497 the α − 1, 4 bonds should be able to take various values. To take this into account, we 498 could randomly pick the dihedral angles of the α − 1, 4 bonds using the Ramachandran 499 plots of their energetically favourable regions. As a first trial, we could use that of 500 maltose, that is well characterised [48]. We expect that, introducing such disorder in 501 the angles, chains will appear longer. Besides, in the model, we consider that all 502 chains are stiff. Thus, to improve the macroscopic representation of the chains, we 503 could introduce the possibility that they bend when encountering steric hindrance. To 504 do so, we would minimise their torsion energy, like it is done in polymer physics 505 models [49]. Opposite to the change suggested for the α − 1, 4 bonds variability, 506 accounting for the flexibility of the chains might lead to a higher granule density, and 507 thus, a lower radius. The fact that these two effects might cancel each other, possibly 508 explains why our simplified model is nonetheless able to capture realistic in vivo 509 granule radii. Besides, in abnormal conditions, the potential formation of double 510 helices may not only lead to glycogen precipitation, but also prevent enzymatic 511 reactions. Thus, a later improvement of the model could include to tune the enzymes' 512 reaction rate depending on the substrate chain configuration and length.

513
Throughout this article, we focus on glycogen synthesis, yet, simulating the 514 degradation dynamics with the algorithms we developed would be straightforward. We 515 expect residual degradation activity to only lightly modify the effective elongation to 516 branching ratio Γ, and slow down the synthesis. In such a case, the CLD would be  In this article, we show that the availability of the substrate strongly influences the 531 enzyme activity, leading to distinct chain lengths and number of non-reducing ends.

532
Noticeably, other modelling approaches, such as metabolic ODE models, instead 533 consider glycogen as a single metabolite approximated by the sum of all glucose units 534 that compose it. A direct consequence is that any structural aspects are neglected and 535 these models cannot differentiate between a single chain of 50,000 glucoses, and an 536 actual granule of the same weight. Still, for instance, the number of non-reducing ends 537 available for elongation are drastically different in these two cases. Thus, coupling our 538 model to glycogen metabolic ODE models, would constitute a hybrid approach that  Different forms of glycogen granules as described by [13]. The first one being aggregates of the second. ρ Radius of a monomer unit described as a sphere in the model. l Distance between two overlapping spheres. Ratio of the elongation and branching reaction rates. R g Radius of gyration (i.e. radius of a glycogen granule in nm) where N is the total number of glucose units in the granule, r i the spatial coordinate of the i th glucose unit, and r mean that of the center of mass of the granule.

A : B
Ratio of the number of A chains (that do not carry any daughter chain) over B chains (that carry at least one). V glucose Volume effectively occupied by all the glucose units incorporated in the granule. V tot Volume of a glycogen granule (i.e. 4/3 · π · R g 3 ) . Ω Occupancy, defined as the ratio of V glucose and V tot . S Score to measure the distance between simulated and experimental CLDs, defined as the mean square difference between the two curves.

Effect of self-exclusion
Modelling the 3D structure of glycogen, while considering self-exclusion of glucoses, allows investigating glycogen chain patterns, molecular density, size of the granules, and crowding at the surface. But is this detailed approach strictly required for studying glycogen's macroscopic properties? In the following, we test various radius values for the spheric description of the glucose units and look whether they affect the macroscopic properties of glycogen. Three scenarios are investigated: ρ = 0 nm, ρ = 0.325 nm, and ρ = 0.650 nm. For ρ = 0 nm, the glucose units have no volume and thus no steric hindrance, allowing chains to overlap. Instead, ρ = 0.325 nm is equal to half of the helix's van der Waals radius, while ρ = 0.650 nm is its total radius. By comparing these scenarios, we can evaluate the impact of steric hindrance on the substrate availability during the synthesis of glycogen. Fig 11 shows that short chains are more abundant for ρ = 0 nm than for the other two values of ρ. Opposite, long chains are more abundant for non-zero ρ. As soon as ρ > 0 nm, since elongation involves adding a single glucose unit, while branching means transferring an entire piece of a branch, it is easier to find the necessary space around the substrate for allowing elongation to take place. This is reinforced by the fact that when adding a new glucose unit at the non-reducing end of a branch, we do not elongate the substrate by the total length of a glucose unit, but only its radial contribution, which is l = 0.24 nm (see Fig 2). As a consequence, steric hindrance stronger impacts branching than elongation. In other words, the number of branching attempts rejected due to steric hindrance is higher than that of elongation. If we would define effective branching and elongation rates that respectively account for these rejections, the branching effective rate would reduce much more than that of elongation. This would lead to an effective elongation to branching ratio Γ eff which increases with the effect of steric hindrance. In the section Elongation to branching ratio, we concluded that as Γ increases, the overall CLD shifts towards higher DPs and the distribution peak decreases. Γ ρ=0 eff < Γ ρ=0.325 eff < Γ ρ=0.650 eff therefore explains both the over-representation of high DPs and the peak reduction as ρ increases in Fig  11. Despite these changes in the CLD, as the steric hindrance increases, we can remark that not only is the number of peak conserved (here unimodality) but their location too.

Flexible location versus strict location branching models
As briefly introduced in the section The model, paragraph Enzymatic reactions, two assumptions for GBE can be made. In the flexible location branching model (used throughout the article) all glucose units that are in the acceptable range can equiprobably be cleaved, and similarly for those that can potentially be transferred. Opposite, in the strict location branching model (introduced here for comparison) GBE always branches at a precise location, i.e. at a given distance from the non-reducing end. For the sake of clarity, these two models and their associated mechanisms are fully detailed in Fig 12. Their patterning impact on the CLD is presented in Fig 13.   Fig 12. Mechanisms of the flexible location and the strict location branching models. Both models fulfill the minimal lengths requirement (Eq 1). In the flexible location branching model, the cleaving position is randomly picked from a uniform distribution, and so is that of branching. Instead, in the strict location branching model, both the cleaving and the branching always occur at a specific distance from the non-reducing end.
In Fig 13, we compare the two branching models, and vary Γ and {L GBE spacing , L GBE leftover , L GBE transferred }, which results in 12 different plots. From the Elongation to branching ratio section, we learnt that increasing Γ spreads the CLD towards higher DPs and reduces the distribution peak typically observed for short chains. In the section Effect of the branching enzyme on the CLD, we observed that increasing the minimal lengths {L GBE spacing , L GBE leftover , L GBE transferred } modifies the shape of the CLD that becomes bi-or even multi-modal. When decreasing these lengths, the multi-modality is weakened while it is reinforced by a small Γ. CLD obtained for small Γ (branching regime), with small minimal lengths for GBE. Top-right: CLD obtained for small Γ (branching regime), with longer minimal lengths for GBE. Increasing the lengths reinforces the multi-modality. Bottom-left: CLD obtained when increasing Γ, with small minimal lengths for GBE. Both conditions contribute to weaken the multi-modality. Bottom-right: CLD when increasing both Γ and minimal lengths for GBE. The first condition tends to weaken multi-modality while the latter one instead reinforces it. As a consequence, multi-modality is observed for the strict location branching model only.
For Γ = 0.1, the strict location branching model (red) matches the flexible one (blue). That is not surprising, since in this particular regime, when branching strongly dominates over elongation, branching occurs as soon as possible (DP > L GBE min ), making the two models equivalent.
When Γ increases (second and third lines), the two models are not anymore equivalent. In the strict location model, the length of the transferred branch is always equal to L GBE transferred , while in the flexible location model, longer chains can be transferred. In both cases, the smallest transferable DP is equal to L GBE transferred (see section 2.2.3). For the strict location model the latter is the most abundant DP, unlike for the flexible one.

Scope of the parameter space
We apply our fitting procedure to all possible combinations of Γ ∈ {0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0, 2.0, 5.0}, L GBE spacing ∈ {1, 2, 3, 4, 5, 6}, L GBE transferred = L GBE leftover ∈ {2, 3, 4, 5}, and ρ ∈ {0 nm, 0.325 nm, 0.650 nm}, while considering two different branching scenarios, i.e. the flexible location and the strict location branching models. It sums up to a total of 1,296 tested parameter sets. The corresponding heat-maps are shown in Fig 14 for the flexible location and Fig 15 for the strict location branching models. In both figures, ρ varies across the columns, while L GBE leftover = L GBE transferred does over the successive rows. On each heatmap, the Y and X axes correspond to L GBE spacing and Γ, respectively. In Fig 14, all scores below 20 are highlighted with white squares, and correspond to good fits (arbitrary cut-off chosen as up to twice the best-fit). Noticeably, various parameter sets fulfill this criterion and all of them show a small L GBE spacing value. As presented in section Comparison to experimental data, paragraph Parameters calibration, the best score is obtained with ρ = 0.65 nm (third column). This supports the fact that steric hindrance plays a role in the chain length distribution of real glycogen, although good matches also exist with ρ = 0 nm, in which the CLD is not impacted by steric hindrance. When focussing on the good scores (dark blue cells), it appears that changing the elongation to branching ratio Γ can be compensated by varying L GBE spacing , or the two other minimal lengths. Opposite, when considering parameter values that are typically reported in the literature, we systematically obtain very poor (i.e. high) scores, even upon varying Γ. This case is highlighted in red and discussed in the section Comparison to experimental data, paragraph Parameters calibration. In general terms, for each heatmap, we observe that L GBE spacing ≥ 4 corresponds to poor scores, i.e. S ≥ 30. This effect is even more pronounced if L GBE spacing ≥ 6, i.e. S ≥ 60 (66.3 and above). Thus, with these parameter sets, our model is not able to reproduce the experimental CLD data by Sullivan and coworkers used for fitting throughout this article [23]. Regarding the enzyme mechanism, this suggests that GBE is able to branch much closer than 4 glucose units away from an existing chain. Similarly, if L GBE transferred = L GBE leftover ≥ 3, we are not able to fit the experimental data by Sullivan and coworkers. Therefore, these minimal lengths must be shorter than those typically reported in the literature.
In Fig 15, we present analogous results for the strict location model. Overall, the fitting procedure returns much poorer scores than for the flexible location branching model (see Fig 14). This is due to the fact that strict branching locations almost systematically lead to a multi-modal distribution, which is not the case of the experimental CLD data set fitted here [23]. Although multi-modality can be cleared by increasing Γ, this would lead to a reduction of the distribution peak and a shift of the overall distribution towards higher DPs. For the experimental data considered, we could not find a good trade-off, corresponding to good fitting scores. As a result, the strict location model seems to be very unlikely, that is why we instead choose the flexible location branching model for this study.

Probabilistic approach to the A:B ratio
In this section, we describe why the A:B ratio must be equal to 1 when no biological nor biophysical properties of the system are considered, but pure probabilities. As presented in Section 2.3.2, we note as two distinct reactions, branching on either an A or a B chain: In this purely mathematical framework, A and B chains have the same chance to react, so that the probability for a reaction to occur only depends on the respective number of A (noted N A ) and B (noted N B ) chains. We sketch the associated probability tree in Fig 16. The tree is symmetric, and the number of paths to a given state follows the binomial coefficients. The horizontal dotted line corresponds to N A − N B = 0, when there are as many A as B chains. We know that in the example of flipping a coin i times, the central limit theorem tells us that the distribution of the difference in the number of heads and tails, tends to a normal distribution centered in 0 when i tends to infinity. This case corresponds to heads and tails having equal probabilities. In our case, a given state of the probability tree (N A , N B ) leads either to the state (N A , N B + 1) with probability p(R A→AB ) = N A N A +N B , or to (N A + 1, N B ) with probability p(R B→AB ) = N B N A +N B . It means that the probability to go from one state to the next one depends on the state of the system, in a way that the probability to come closer to N A − N B = 0 (horizontal dotted line) is always higher than that of spreading away. Additionally, the distribution of the probabilities remains symmetric with respect to the case N A = N B . Based on these considerations, the mean value N A − N B = 0 and the distribution is even more peaked than in the simple case of flipping a coin. Since the mean value N A − N B = 0, the mean ratio N A N B = 1, in other words the A:B ratio is 1.