**Authors:**

Yvan Rousset, et al.

**Citation: **Rousset Y, Ebenhöh O, Raguin A (2023) Stochastic modelling of a three-dimensional glycogen granule synthesis and impact of the branching enzyme. PLoS Comput Biol 19(5):

e1010694.

https://doi.org/10.1371/journal.pcbi.1010694

**Editor: **Mark Alber,

University of California Riverside, UNITED STATES

**Received: **October 28, 2022; **Accepted: **March 25, 2023; **Published: ** May 19, 2023

**Copyright: ** © 2023 Rousset et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability: **The source code and data used to produce the results and analyses presented in this manuscript are available on a Gitlab repository at https://gitlab.com/qtb-hhu/models/Stochastic-modeling-of-a-three-dimensional-glycogen-granule. The version associated to this publication is tagged as ”v1-updated”.

**Funding: **This paper is supported by European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement PoLiMeR, No 812616, that funds the position of YR. OE and AR are supported by the Deutsche Forschungsgemeinschaft (DFG) under Germany’s Excellence Strategy EXC 2048/1, Project ID: 390686111. The current position of AR is funded by the German federal and state programme Professorinnenprogramm III for female scientists. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests: ** The authors declare that no competing interests exist.

## 1 Introduction

Management of energy resources is crucial for an organism to survive, since nutrient availability may vary considerably with time. Moreover, organisms face numerous other stresses that may temporarily increase energy demand. To react to such changes in energy supply and demand, it is apparent that some internal stores are necessary. These stores are filled when nutrients are abundant and depleted when demand exceeds available supply. As direct substrate for catabolic pathways, glucose plays a central role in energy metabolism in most organisms [1]. While plants have evolved to store glucose as starch, animals, fungi, and most bacteria store glucose as glycogen. Both starch and glycogen are branched polymers consisting of glucose monomers, linked into linear chains by *α*-1,4 bonds, and branching points by *α*-1,6 bonds. However, these two macromolecules exhibit rather different physico-chemical properties. In contrast to glycogen, starch is insoluble in water under physiological conditions, and contains high density crystalline regions. These different properties are reflected by distinct branching patterns and chain length distributions (CLD). Functionally, starch and glycogen can be compared to capacitors in electric circuits. The latter are able to store and release electrons depending on current and voltage. Thus, they can be used to stabilise a fluctuating electric signal. Analogously, glycogen and starch can be seen as two different capacitors that both contribute to glucose homeostasis by managing energy resources through time.

While the fine structure of starch has been widely investigated over the past 50 years, less is known on that of glycogen [2]. The characterisation of the detailed structure of glycogen, as well as the interplay between its structural properties and the dynamics of glycogenesis (synthesis) and glycogenolysis (degradation) is unclear. Yet, both glycogen structure and dynamics are of utmost interest for understanding glycogen metabolism and the impact of related genetic variances. For human health, this is particularly important considering the increasing prevalence of glycogen storage diseases (GSDs), as well as diabetes, and other glycogen related disorders.

So far, different structures of glycogen have been proposed [3–6], but a structural model that emerges from the detailed underlying enzymatic mechanisms of synthesis and degradation is still lacking. Understanding and precisely describing the biochemistry of glycogen is challenging. With a molecular weight of 10^{6} to 10^{7} g⋅mol^{−1} [7–9], glycogen is a large molecule, even compared to the enzymes involved in its dynamics [10–12]. It is thus expected that, enzymes synthesising, degrading, or otherwise altering glycogen, can only access certain branches near its surface, while many glucose residues near the centre of the molecule remain ‘hidden’. At the macroscopic scale, the structure of glycogen is well known. Drochmans [13] observed two populations of glycogen granules in rat liver. The so-called *α* granules are aggregates of the smaller *β* granules [14–16]. The latter have a radius in the range from 10 to 20 nm, while the radius of *α* granules ranges between 40 and 300 nm [13]. Here, we focus on *β* granules, whose synthesis and degradation both involve a relatively small number of enzymes (see Fig 1). For a comprehensive review, see [17].

Fig 1. Main enzyme reactions involved in the synthesis and breakdown of glycogen.

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

To initiate a new glycogen molecule, a chain of 5 glucose units is synthesised and bound to a glycogenin (GN) protein, which forms the core of the final granule [18–21]. Once initiation is completed, the granule is expanded by the two enzymes Glycogen Synthase (GS) and Glycogen Branching Enzyme (GBE). GS is an elongating enzyme that adds one glucose residue to the non-reducing end of an *α*-1,4 linear chain, thereby forming a new *α*-1,4 glucosydic bond. GBE cleaves a small part of a linear chain and creates a new branch by forming an *α*-1,6 glucosydic bond. We call “daughter” the newly formed chain and “mother” the one it is branching from.

Besides synthesis, granules are subject to degradation, that is performed by two other enzymes. Glycogen Phosphorylase (GP) and Glycogen Debranching Enzyme (GDE) respectively shorten and debranch glycogen branches.

Depending on the relative kinetics of these four enzymes (GS, GBE, GP, and GDE) the overall size of the glycogen granule can either increase or decrease. In this article, we choose to focus on glycogen synthesis, and more specifically on the role of the branching enzyme GBE.

Experimental observations [22, 23] of glycogen show an average chain length (CL) around 13 glucose units, which depends on the organism type. A typical peak is observed at low degree of polymerisation (DP), around DP 8, and almost no chains are detected above DP 40 [22–25]. The degree of branching is defined in two ways in the literature. It is most commonly defined as the ratio of *α*–1,6 to *α*–1,4 linkages, but sometimes also as the average number of *α*–1,6 bonds per chain [6]. We will apply the first definition throughout the paper. This ratio is in the range 0.02 − 0.05 in amylopectin (the prime constituent of starch) [26–28] and 0.06 − 0.10 in glycogen [9, 29]. In 1956, Peat et al. [30] introduced the concept of A and B chains. An A chain does not carry any branch, while a B chain does. The A:B ratio is an important characteristics of glycogen granules and an indicator of the branching pattern. Early studies reported an A:B ratio of 1 in glycogen [22, 31], while it is usually greater than 1 in amylose, and ranging from 1.5 to 2.6 in amylopectin [31, 32].

As early as in the 1940s, various hypotheses have been formulated aiming at explaining macroscopic features of glycogen granules. One of them has become known as the Whelan model, which assumes that every chain carries on average two branches. Based on the Whelan model, the nowadays widespread idea emerged that glycogen can be described as a fractal molecule [6, 33, 34]. A fractal glycogen structure is indeed attractive because it can reproduce various structural properties of glycogen. Moreover, it provides a mechanism explaining why glycogen granules seem to have a maximum size. With this model, glycogen becomes too dense at the surface due to an exponential increase of the number of non-reducing ends with the distance from the centre. Thus, steric hindrance prevents synthesis to continue. However, Manners [2] stressed in 1990 that there is no clear evidence that supports a regular branching and therefore a fractal pattern. More recently, further arguments and results against a fractal structure have been raised [35–37]. Besides, independently from any fractal structure considerations, Zhang et al. [38] proposed a mathematical model based on a Monte Carlo approach to numerically simulate glycogen biosynthesis, aiming to support the common idea that steric hindrance limits granule growth. In this model, glucose units are placed on a three-dimensional grid, and the granule biosynthesis is simulated by adding glucose units on any of the 26 neighbouring positions around the end of the growing chain. As a result, limited growth is an emergent property of this model.

Here, we propose a mechanistic model for glycogen synthesis, focussing on the impact of the branching enzyme on the granule structure. We aim at explaining macroscopic and experimentally observable quantities, such as chain length distributions (CLDs), from the underlying enzymatic mechanisms. Such CLDs can be predicted theoretically using numerical inversion of Laplace transforms [39] or kinetic equations [40–42]. However, besides the fact that CLDs are only one of the many quantities that can be measured in the glycogen structures we simulate, our model accounts for complex features for which an analytic treatment is no longer feasible. For instance, these are the complete connectivity of the structure, complex enzyme mechanisms, the helical structure of the linear chains, and steric-hindrance effects. Also, in contrast to the model proposed by Zhang et al. [38], we do not restrict the location of glucose units to a grid, and instead reflect the helical structure of glucose chains, using parameter values derived from biophysical properties. Therefore, our model provides a more flexible and more realistic representation of the three-dimensional granule structure. Specifically, it is designed to study the effect of the enzymatic activity on the structure, and infer model’s parameter values using comparisons to quantitative experimental data. Differently, Zhang et al. [38] focus on explaining the limited growth of glycogen granules. To our knowledge, these scientific questions have so far not been addressed by computational approaches. In this article, we first detail the geometrical and structural features, and enzymatic reactions, taken into consideration in our model. Then, we analyse distinct properties of the model with a specific focus on the effect of the branching enzyme. Finally, we compare the model to experimental data and discuss the parameter values resulting in a best fit, in relation to typical values reported in the experimental literature. In addition, several complementary results justifying our modelling choices are reported in the extensive Supporting information (S1 Text–S5 Text).

## 3 Material and methods

### 3.1 Simulation procedure

We present a model that simulates the dynamics of glycogen synthesis. We record all enzymatic reaction steps, the time at which they occur, and the full three-dimensional details of the glycogen structure at any time point. The model specifically keeps track of each glucose position, the complete connectivity of the chains, and the position of each branching point. To account for this complexity, we implement the model using a stochastic algorithm. This approach also allows to specifically consider how changes in the glycogen structure enable or disable enzymatic reactions. For comparison, employing a deterministic approach, for instance based on systems of ordinary differential equations, would lead to unnecessarily complex simulation rules, that would also include additional ad-hoc assumptions. Besides, randomness has a stronger impact on a system as it involves small numbers, like the one we consider here. Indeed, the synthesis and breakdown of a single glycogen granule of ca. 50,000 glucose units involves only a small number of enzymes. For example, experiments indicate that on average a single glycogen synthase enzyme is active per granule [53]. As a result, a stochastic approach appears very natural to deal with this complex and spatially structured substrate.

We base our simulations on the Gillespie algorithm [54, 55]. At each time step, it consists in randomly selecting both an enzymatic reaction and its duration, while systematically updating the structure of the glycogen molecule accordingly. Although the enzymatic reaction is selected from a uniform distribution, the time is chosen from an exponential distribution, such that the Gillespie algorithm allows to simulate the real time of the dynamics of the system, as far as the underlying enzymatic mechanisms and their kinetic parameter values are known. Each reaction is associated to a specific propensity, that we choose as proportional to the concentration of enzymes and substrate, thereby following a typical mass-action kinetics.

The Gillespie algorithm accounts for all possible reactions, and keeps track of any modification of the granule structure, such that only possible reactions can be selected. Nonetheless, because of the complexity of the structure we simulate, to account for steric hindrance among different branches, we instead allow for certain reactions to be rejected. For these rejected reactions, the time elapsed is not accounted.

We provide details about the Gillespie algorithm (Fig A in S5 Text), and how it was employed to simulate the dynamic changes of the three-dimensional structure of the granule (Fig B in S5 Text). Moreover, the source code of our model, together with Jupyter Notebooks that recreate the main figures of this manuscript, are available on our gitlab repository (link provided below under Data availability).

### 3.2 Best fit algorithm

The model contains various kinds of parameters. Some describe the physical properties of the glycogen structure, others relate to the enzymatic activity, including the enzyme substrate specificities. On the one hand, certain parameter values are inferred from literature data, for instance the minimum DP for GS to act (). On the other hand, other parameter values are free to be fitted by our custom designed best fit algorithm. Here we choose to fit the minimum DP between two branches after a branching reaction (), the minimum DP between the non-reducing end of the mother branch and the new branching point (), and the ratio between the elongation and branching reaction rates (Γ).

Our best fit algorithm consists in setting bounding ranges for the values of the parameters to be fitted, and systematically scanning the parameter space thereby defined. For each set of parameter values tested, we run 50 simulations, take the average of the resulting CLDs, and compare it to our targeted experimental data. For each set of parameter values, the comparison to the experimental data is measured by the score , defined as follows:

where is the experimental CLD abundance for a given DP and the average simulated CLD abundance for a given DP. Therefore, measures the squared difference between the average simulated and the experimental CLDs. The best fit is found for the set of parameter values that has the lowest score .

## 4 Conclusion

### 4.1 Discussion of the results

According to the World Health Organization, metabolic diseases are a burden in western countries. Among them, Glycogen Storage Diseases (GSDs), Lafora disease, Adult Polyglucosan Body Disease (APBD), and even diabetes directly or indirectly involve glycogen. Investigating the regulation of glycogen’s structural properties could therefore strongly contribute to further understand such diseases. Computational models that encompass the complex interplay of both glycogen’s structure and its metabolism, allow to tackle this challenge but remain poorly exploited.

In this article, we present a stochastic structural model for glycogen synthesis. The model provides a precise description of both the structure of glycogen in three dimensions, that we can visualise, and the detailed dynamics and mechanistics of the underlying enzymatic reactions. For instance, both elongation and branching have precise rules regarding their substrate specificities. Modelling glycogen granules in three dimensions is made possible by a coarse-grained geometrical description at the glucose level, allowing us to track all the glucose units in space. This description also accounts for the steric hindrance effects resulting from the impossibility for the chains to overlap. Simulating reaction dynamics relies on a Gillespie algorithm, which determines when and where a reaction of branching or elongation takes place, thereby dictating the corresponding change in the three-dimensional structure.

We show qualitatively how enzyme activity affects glycogen structure for generic sets of parameter values. We highlight two different synthesis regimes, depending on Γ, the ratio between the elongation and branching reaction rates. By varying this ratio, either small and dense, or big and sparse granules can be simulated. *In vivo*, it can be expected that Γ depends on numerous factors, such as the organism under investigation, the cell type, and possibly the external conditions. In addition, a model that would consider chains bending and intermolecular interactions might lead to more complex results, in particular during the synthesis of sparse granules. In our results, the phenomenology of the two synthesis regimes is also confirmed by the occupancy profiles along the radius of the corresponding granules. In both synthesis regimes, first the center of the granule is filled around GN, before reaching a critical density, preventing further internal reactions to occur due to steric hindrance. Then, the granule expands such that the density remains approximately constant within the sphere defined by the gyration radius. While this result is a consequence of the geometrical assumptions of our model, it is interesting to note that a glycogen fractal description would instead give a density exponentially increasing with the granule radius. Besides our own results and the various arguments exhibited against a fractal representation of glycogen, we would like to highlight here that, as soon as one considers even a single degree of freedom in the glycogen branching reaction, it would rapidly lead to the loss of any “fractal-like” structure. Such degrees of freedom are necessarily present in natural conditions. For instance, the dihedral angles defining the *α* − 1, 6 bonds may take various values, making a fractal pattern very unlikely *in vivo*. Beyond these spatial considerations, the model establishes a natural and clear connection between enzymes’ reaction rates, and both the degree of polymerisation of the chains, and the number of non-reducing ends.

Our results show that the chain length distribution of glycogen is highly sensitive to the branching reaction, predominantly its mechanism. Each of the three characteristic minimal lengths of the reaction has specific effects on the CLD. Additionally, if any of them increases, less branching outcomes are possible, eventually leading to a bi-modal or even multi-modal distribution. This effect is enhanced by high branching activity. When varied together, these minimal lengths show even more complex imprints on the CLD. In contrast, multi-modal distributions become less pronounced by increasing the elongation reaction rate, because rapid elongation increases the number of possible configuration outcomes. Additionally, increasing elongation leads to longer chains and results in a CLD spreading towards higher DPs. Altogether, we show that the CLD, and in particular peaks’ location and intensity, are subtly affected by several complex effects.

Guided by these findings, we propose to consider the CLD not only as an important structural feature of glycogen, but also as a signature of GBE. Thus, fitting experimental data with our model arises as a natural strategy to infer knowledge on the GBE mechanism. Not only did we illustrate the strength of our fitting procedure on an experimental data set from Sullivan et al. [25], but we also extracted several macroscopic characteristic features from the resulting fit, which we compared to various literature sources, thereby confirming the validity of our results. Using this fit, at the microscopic level, we were able to discriminate between the two branching models we hypothesised, and selected the flexible location branching over the strict one. Besides, we could critically evaluate parameter values typically reported in the literature. For instance, it is often assumed that GBE transfers branches of DP 4 or 6, or even longer [33, 49, 50]. Similarly, it is typically reported that around 6 glucose units space two branches. Within the framework of our model and its underlying assumptions, it is impossible to reproduce CLDs of *in vivo* glycogen using the above mentioned values. Instead, our fitting procedure suggests that a high flexibility is necessary for both the branching mechanism and the minimal lengths involved. This finding fully confirms the importance of modelling glycogen synthesis using a stochastic approach.

Moreover, our customisable branching model shows that the A:B ratio is independent of the kinetic parameter Γ, but specifically determined by the difference between and . Based on these observations, we hypothesise that the branching mechanism is chiefly responsible for the structural differences observed between starch and glycogen.

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

### 4.2 Outlook

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

Throughout this article, we focus on glycogen synthesis, yet, simulating the degradation dynamics with the algorithms we developed would be straightforward. We expect residual degradation activity to only lightly modify the effective elongation to branching ratio Γ, and slow down the synthesis. In such a case, the CLD would be slightly shifted to the left, and the first chain length detected would be the smallest value between and its GP counterpart. It could be particularly interesting to account for the degrading enzymes, beyond a residual activity, in case the synthesis is defective. Indeed, abnormal structures produced by defective synthesis enzymes could be degraded, and thereby corrected, by degrading enzymes. For instance, if becomes too short for GS to elongate the newly formed chain, GDE could unbranch the latter and thereby preserve the macroscopic properties of the granule. Last, following an approach analogous to the one taken in this article, one could choose to investigate the glycogen granules’ breakdown in full depth, by first synthesising granules and then proceeding to their degradation. Although, for the sake of simplicity, we would then uncouple in time the synthesis and degradation processes, it would still be very interesting to study the mutual impact of distinct modes of synthesis and degradation on the overall glucose release and fixation.

In this article, we show that the enzyme substrate specificity strongly influences the enzyme activity, leading to distinct chain length distributions and number of non-reducing ends. Noticeably, other modelling approaches, such as kinetic models using systems of ordinary differential equations (ODEs), instead consider glycogen as a single metabolite, approximated by the sum of all glucose units that compose it. A direct consequence is that any structural aspects are neglected and these models cannot differentiate between a single chain of 50,000 glucoses, and an actual granule of the same weight. Still, for instance, the number of non-reducing ends available for elongation are drastically different in these two cases. Thus, coupling our model to glycogen metabolic ODE models, would constitute a hybrid approach that would include key structural details, while enlarging its biochemical scope. It would thereby open up a whole new range of modelling possibilities. For instance, it would allow to investigate the interplay between glycogen structure and the evolution in time of important metabolites (e.g. glucose-1-phosphate, glucose-6-phosphate, and UDP-glucose), under various physiological conditions, including diseases scenarios. Using this approach, we shall be able to characterise the phenomenology of each glycogen storage disease, with a focus on the role of glycogen structure, and address questions such as glycogen accumulation, glucose cycling, glucose homeostasis, and even glycogen precipitation.