Unravelling heparin's enhancement of amyloid aggregation in a model peptide system

Abstract

A coarse-grained (CG) model for heparin, an anionic polysaccharide, was developed to investigate the mechanisms of heparin's enhancement of fibrillation in many amyloidogenic peptides. CG molecular dynamics simulations revealed that heparin, by forming contacts with the model amyloidogenic peptide, amyloid-β's K₁₆LVFFAE₂₂ fragment (Aβ_16–22), promoted long-lived and highly beta-sheet-like domains in the peptide oligomers. Concomitantly, heparin-Aβ_16–22 contacts suppressed the entropy of mixing of the oligomers’ beta-domains. Such oligomers could make better seeds for fibrillation, potentially contributing to heparin's fibril-enhancing behaviour. Additionally, reductions in heparin's flexibility led to delayed aggregation, and less ordered Aβ_16–22 oligomers, thus offering insights into the contrasting inhibition of fibrillation by the relatively rigid polysaccharide, chitosan.

---

1 Introduction

Glycosaminoglycans (GAGs) are linear anionic polysaccharides found abundantly in physiological environments like the extracellular matrix.^1,2 In addition to sharing these microenvironments with amyloid-forming peptides, GAGs are well known to co-aggregate with amyloid fibrils.^3,4 Many GAGs, including the highly anionic model GAG heparin, are also known to enhance the fibrillation of peptides including the Alzheimer's disease-related amyloid-β (Aβ),^5,6 several disease-related peptides,^7–9 and peptide hormones.^10–13 Beyond merely enhancing fibrillation, heparin has also been shown to induce specific fibril polymorphs and trigger condensate-to-fibril transitions.^14,15 Yet, the mechanisms by which heparin and other GAGs modulate peptide aggregation remain elusive.

An interesting contrast to the enhancement of amyloid aggregation by heparin is the inhibition of amyloid aggregation by another pyranosic linear polysaccharide, chitosan. Where heparin is flexible and polyanionic, chitosan is relatively rigid and displays pH-dependent polycationicity. These opposite effects with heparin and chitosan have been reported for several peptides including Aβ₄₂ (net charge −3 under physiological conditions),^16,17 tau_298–317 (net charge +3 under physiological conditions),⁷ α-synuclein (net charge −9 under physiological conditions)⁸ and amylin (net charge +6 under physiological conditions).¹⁸ The diverse net charges displayed by these peptides, from −9 to +6, show that the opposite effects of heparin and chitosan on their aggregation occur independently of the peptides’ net charge. However, the mechanistic details of this phenomenon are not well characterized.

Atomistic molecular simulations, while adept at tackling such gaps in mechanistic understandings, have so far been limited to small numbers of peptides^19–21 and short timescales.²² This is a pervasive limitation of classical atomistic forcefields, where rugged free energy landscapes lead to poor sampling and high computational costs. Much of these costs can be overcome with coarse-grained (CG) molecular dynamics simulations, whereby groups of atoms are coarsely represented by their average chemical characteristics as “beads”.

The forcefield for CG molecular dynamics developed in our lab, ProMPT, allows for the unbiased folding of secondary and super-secondary protein structures while retaining an explicit solvent architecture.²³ In this work, we developed a ProMPT parameter set for the model GAG heparin (Fig. 1(a)), which was found to closely match the characteristic torsional angles of heparin in atomistic simulations, and experimentally measured radii of gyration.²⁴ The development of this heparin model thus enabled us to study the effects of GAGs on amyloid aggregation using the ProMPT forcefield for CG molecular dynamics.

Fig. 1 (a) Mapping of the CG heparin model. (b) Schematic representation of the CG Aβ_16–22 fragment K₁₆LVFFAE₂₂ (left). A key for the bead types, and a schematic depiction of the CG MARTINI polarizable solvent model (right). The anionic and cationic dummy particles of any polar bead X would be named Xm and Xp, respectively.

We chose the extensively studied amyloidogenic fragment of Aβ, K₁₆LVFFAE₂₂ (referred to henceforth as Aβ_16–22),²⁵ as the model peptide for this study. The Aβ_16–22 fragment spans the heparin-binding basic patch of full-length Aβ at K₁₆,^5,26 thus making it suitable for our goal of studying the role of heparin in amyloid aggregation. Fig. 1(b) shows a schematic representation of Aβ_16–22 as defined in ProMPT.

Through CG molecular dynamics simulations, we studied how the presence of heparin, heparin's degree of polymerization (dp), and rigidity influenced the aggregation of ordered Aβ_16–22 oligomers. We also discussed our results on the effects of heparin's rigidity on Aβ_{16– 22}'s ordered aggregation in the context of previous studies, including one from our group,¹⁷ on chitosan's inhibition of Aβ_16–22 aggregation.^5,16

2 Methods

2.1 CG forcefield details

Protein model. Aβ_16–22 peptides were modelled with the ProMPT forcefield.²³ (Fig. 1(b) left). Bonded and non-bonded parameters for all amino acids have been detailed in the corresponding work of Sahoo et al.²³ While ProMPT allows the implementation of Cα-style backbone dihedral potentials for secondary structure biasing, no such biasing potentials were applied in this work, in step with our previous studies with Aβ_16–22.¹⁷

Heparin model. Heparin was modelled by repeats of the most common heparin disaccharide, [α-L-iduronic acid-(1 → 4)-α-D-glucosamine-(2,6)-disulfate], abbreviated here as (IDO-GDS). The atoms were grouped into CG beads such that their partial charges would sum up to either 0 e or −1 e. In effect, the structure of each CG monosaccharide is represented by four main beads: anionic beads B2 and B6 representing sulfate and carboxylate groups on the 2′ and 6′ carbon positions of the pyranose ring, and polar beads BG and B3 representing the glycosidic backbone, and 3′ carbon group atoms with net-zero charge (Fig. 1(a)). The length of a heparin chain is defined by its degree of polymerization (dp), e.g., a heparin chain of dp18 is composed of eighteen monosaccharides or nine (IDO-GDS) repeats.

While heparin is a polydisperse molecule (with molecular weights of up to 30 000 Da depending on the source²⁷), low molecular weight heparin (LMWH) is a formulation with a consistent average molecular weight of ∼5000 Da, i.e., heparin chains roughly dp18 in length.²⁸ Therefore, most of our simulations employ dp18 heparin—a stand-in for laboratory/pharmaceutical grade LMWHs.

Further details on the parametrization of the model's nonbonded and bonded interactions, validation of the model against experimental values of radii of gyration, and characteristic torsion data from atomistic simulation are presented in Section 1 of the ESI.†

Solvent model. The MARTINI polarizable water model²⁹ was used as the solvent in our simulations. Standard MARTINI monovalent ions were used to balance net charge where necessary. Schematic structures are presented in Fig. 1(b) (right).

A brief note on the salient differences between ProMPT and MARTINI. The ProMPT forcefield is based on the MARTINI 2.2P parameter set for biomolecules,³⁰ but has some important distinctions that enable the folding of proteins into secondary and super-secondary structures. The first difference is the introduction of structural polarization in the form of internal Drude-like oscillating dummy charges, which emulate the inherent dipole moments of peptide bonds, and polar amino acid side chains. All polar beads in heparin and Aβ_16–22, indicated by pink and purple in Fig. 1, contain such internal oscillating charged dummies. In a way similar to the MARTINI polarizable water model,²⁹ the Lennard-Jones potentials of the polarized beads are scaled down to correct for the additional electrostatic interactions via the internal dummy charges. The second difference is the enhanced hydrophobicity of the side chains of the apolar amino acids. Previous publications contain details on the parameterization of ProMPT, its predecessor forcefield WEPPROM, and studies on how these alterations assist in the folding of proteins, with optional secondary-structure specific Cα-dihedral potentials in the case of ProMPT and without any Cα-dihedral bias in the case of WEPPROM.^23,31,32 For a complete summary of Lennard-Jones interactions for the bead types used in this work, see Table S2 (ESI†).

2.2 CG simulation details

Starting with a fully extended conformation, a single Aβ_16–22 peptide was placed in a water box and subjected to 10 000 steps of steepest descent energy minimization and 10 ns of NPT simulation to generate a random initial structure. The coordinates of Aβ_16–22 from the final frame, a collapsed conformation without helical or beta-sheet-like secondary structure, served as the initial configuration when building the simulation cell. Similarly, starting with the CG-mapped structure of the PDB record 3IRI,³³ a single dp18 heparin was placed in a water box and subjected to 10 000 steps of steepest descent energy minimization and 10 ns of NPT simulation to generate a random initial structure (the dp18 heparin structure was clipped to obtain coordinates for shorter heparin fragments). These coordinates of Aβ_16–22 molecules and heparin molecules were randomly inserted, at least 1 nm apart, in a 9 × 9 × 9 nm cubic periodic box to initialize each system. Each system was solvated with roughly, 7600 particles of MARTINI polarizable solvent, and monovalent ions to balance net charges where necessary.

These initial system configurations were equilibrated with 10 000 steps of steepest descent energy minimization, followed by 50 000 steps of NPT simulation at 0.01 ps timesteps keeping the positions of solute molecules restrained with spring potentials. Finally, 3000 ns of NPT production MD was performed with timesteps of 0.01 ps. Four independent trials were performed, each 3 μs long and equilibrated with a unique velocity seed.

We used the leapfrog integrator in conjunction with the Nose–Hoover thermostat at 350 K with a time constant of 1 ps.³⁴ Solvent, heparin, and Aβ_16–22 molecules were coupled to separate temperature baths. Pressure was maintained at 1 bar with an isotropic Parrinello–Rahman barostat, 5 ps time constants, and compressibility of 3.5 × 10⁻⁵ bar⁻¹.³⁵ Long-range electrostatics were computed with the Particle Mesh Ewald scheme with a relative electrostatic permittivity of 2.5.³⁶ Neighbour lists for short-range interaction calculations were updated every 10 steps. LINCS was used to constrain the dummy bonds within the MARTINI polarizable solvent particles, and the bonds within the aromatic rings of the protein molecules.³⁷ All simulations were performed using the GROMACS 2019.4 simulation engine.³⁸

2.3 Analysis

Definitions of Aβ_16–22 peptides, aggregates and consolidated aggregates. An Aβ_16–22 aggregate is defined by two or more peptides with at least one inter-peptide contact between their non-dummy particles (0.7 nm cutoff). An aggregate was considered “consolidated” if it contained all the peptides present in the simulation box. The term “peptide” is used here to refer to any Aβ_16–22 molecule at any state—monomeric or otherwise.

Definition of beta-domain. A “beta-domain” is a region of local order in an aggregate, defined as a series of at least three peptides, each connected by a stretch of 4 or more contiguous BB–BB contact pairs. The peptides constituting a beta-domain are called “beta-strands”. Beta-domains are identified by the sets of their constituent peptide id's, e.g., {1,3,4,8}, and {2,5,9} could constitute distinct beta-domains.

Beta-domain lifetime. Lifetimes of beta-domains were calculated by counting the number of nanoseconds, generally in the last 1 μs of simulation, for which a beta-domain (identified by its constituent peptide id's) was observed. Lifetimes were averaged by weighing lifetimes by themselves so that long-lived beta-domains are given more importance than very short-lived domains with near-zero lifetimes.

Order parameter, Q. End-to-end vectors of beta-strands, defined between the BB beads of K₁₆ and E₂₂, were used to compute the orientational order parameter, Q, for each beta-domain. The order parameter Q is given bywhere θ is the angle between a peptide's end-to-end vector and the director vector indicating the preferred local direction.

Heparin contacts per beta-strand H_beta. A contact between a beta-strand and heparin was defined by a distance of 0.7 nm or less between BB, S1, S2 or S3 beads of the beta-strand and BG, B2, B3 or B6 beads of heparin. The average number of heparin contacts per beta-strand, H_beta was determined by counting the number of heparin contacts of a beta-domain, averaged over its lifetime and divided by the number of beta-strands in the domain.

Number of beta-strands, N_beta. N _beta refers to the number of beta strands at any given instant of time. It is calculated by summing over the number of beta-strands across all the beta-domains present at any given time. N_beta is reported as an average across four replicas.

Number of beta-domain combinations, N_comb. N _comb is the sum of unique beta-domain combinations, defined by the set of the peptide id's, observed over a period of time in a single simulation. N_comb is reported as an average across four replicas.

Compactness of Aβ_16–22 aggregates. Compactness refers to the ratio of the smallest and largest moments of inertia of an aggregate. The smallest and largest moments were extracted from the moment of inertia tensor, calculated using the moment of inertia method of the AtomGroup class of the MDAnalysis python package.³⁹

Characteristic aggregation time. To characterize the time taken for the consolidation of all peptides into a single aggregated cluster, we constructed a binary time series given by g(t):

A characteristic aggregation time, τ_agg was calculated from this time series by fitting a sigmoid curve:

3 Results and discussion

First, we investigated how a single heparin chain, with a degree of polymerization of 18 (dp18), would influence Aβ_16–22 aggregation at different concentrations—16, 24 and 32 mM, mapping to numbers of peptides (N_pep), 8, 12 and 16, respectively in a 9 × 9 × 9 nm periodic solvent box.

In each simulation, a single consolidated aggregate comprising all the available Aβ_16–22 peptides was obtained. Representative structures of consolidated aggregates in water and with heparin at N_pep = 16 are shown in Fig. 2(a) and (b), respectively. Aβ_16–22 oligomers with characteristic hydrophobic cores were obtained in both cases (Fig. S7, ESI†). Conforming to its spheroid shape, heparin bound to the periphery of the Aβ_16–22 oligomer, primarily with the cationic side chains of Aβ_16–22 's K₁₆ residue, and secondarily with peptide backbones though the B6 and B2 beads of heparin's GDS and IDO subunits, respectively (Fig. S8, ESI†). The co-existence of electrostatic interactions between heparin and K residues, and polar interactions of heparin with the backbones of diverse amino acids was also reported in a recent atomistic simulation study of heparin with the R3 fragment of the tau peptide.²¹ It is well-known from past experimental studies that heparin chains are structurally integral to mature amyloid fibrils.^5,6,14 Our results suggest that heparin may also be an integral component of pre-fibrillar oligomers.

Fig. 2 (a) and (b) Representative snapshots of Aβ_16–22 aggregates in water (a) and with dp18 heparin (b) at N_pep = 16. Beta-domains are coloured in red and blue, while all other peptides are coloured in pink, and heparin is coloured in orange. (c)–(e) Average number of beta-strands (N_beta), violin plots of distributions of Q, and average number of beta-domain combinations, N_comb, respectively, for Aβ_16–22 in water (pink, ■ dotted lines) and with dp18 heparin (orange, •, solid lines) at varying N_pep. (f)–(h) Probability distributions of beta-domains’ Q against their lifetimes for Aβ_16–22 aggregation at N_pep = 16 in water (f) and with dp18 heparin (g). Probability distribution of heparin contacts per peptide, H_beta, and lifetimes of beta-domains for Aβ_16–22 at N_pep = 16 with dp18 heparin (h).

While we cannot access the timescales and system sizes required to sample fibrillar structures, we can study the fibril-like qualities of early oligomeric species. A key feature of amyloid fibrils is an enhancement in beta-sheet secondary structures. In our simulations, Aβ_16–22 oligomers contain transient patches of beta-sheet-like local order characterized by four or more contiguous backbone–backbone contacts between three or more peptides (such as the red and blue peptide domains in Fig. 2(a) and (b)). We call such patches “beta-domains”, and their constituent peptides “beta-strands”. While adept at detecting locally aligned peptides within oligomers, our beta-domain concept cannot discriminate between collapsed and extended structures. Consequently, the average number of beta-strands, N_beta, in water and with heparin were nearly identical (Fig. 2(c)). Unless specified otherwise, data for N_beta and all metrics in all figures were reported from the final 1 μs of four independent replica simulations, and averages were reported with error bars of ±2 standard errors.

Looking beyond raw beta-strand counts, we looked for more subtle structural differences in the orientation of beta-strands within the beta-domains with the order parameter, Q. A Q value of 1 corresponds to a set of peptides oriented in a perfectly parallel (or antiparallel) manner, while a value of 0 corresponds to randomly oriented peptides. In practice, highly ordered beta-sheets display Q values in the 0.75–0.95 range.⁴⁰

Relative to Aβ_16–22 in water, Fig. 2(d) shows a consistent increase in populations of beta domains with Q ≥ 0.75 in the presence of a single chain of dp18 heparin across all peptide concentrations studied. Therefore, we concluded that heparin enhanced fibril-like order in oligomers, a trend that agrees with the enhanced fibrillation reported in previous experimental studies.^5,41

We could also track how the different peptide strands, identified by their peptide id in the simulations’ topologies, combined and recombined to form beta-domains. An example of two beta-domains recombining would be the domains identified by peptide id's {1,3,4,8}, and {2,5,9} recombining into two new domains {1,3,4}, and {2,5,9,8}. Importantly, by our definition, {1,3,4,8} and {8,1,4,3} are identical combinations, i.e., the order of peptide id's is irrelevant.

Fig. S10 (ESI†) illustrates, over a 300 ns period of self-assembly, the differences in the propensities for beta-domain recombination—rampant in Aβ_16–22 aggregation in water, and significantly retarded in Aβ_16–22 aggregation in with dp18 heparin. These differences were quantified by the numbers of beta-domain combinations, N_comb, and their lifetimes. Across different N_pep, N_comb in water far exceeded that in the presence of heparin (Fig. 2(e)). Complementarily, the mean lifetimes of the beta-domains were higher in the presence of heparin than without (Fig. S11, ESI†). These data indicated that the entropy of mixing of beta-domains, S^mix = k_B ln(N_comb) where k_B is the Boltzmann constant, was suppressed by heparin.

Differences in the beta-domains’ mixing propensities could be traced back to their level of order and their interactions with heparin. At N_pep = 16, beta-domains with long lifetimes generally had higher Q (Fig. 2(f) and (g)), but heparin significantly shifted the ensemble of beta-domains towards Q ≥ 0.75, and lifetimes over 100 ns (see ESI† for data at other N_pep values). Long-lived beta-domains also had, on average, more heparin contacts per beta-strand (H_beta), indicating that beta-domain mixing was suppressed by the energetic costs of breaking peptide–heparin contacts (Fig. 2(h)). The following mechanism was summarized: peptide–heparin contacts prolonged beta-domain lifetimes, which then allowed the constituent beta-strands to arrange themselves into highly ordered configurations. By this mechanism, heparin could potentially spawn highly ordered oligomeric seeds for fibrillation, which could explain heparin's fibril-enhancing properties for Aβ and other peptides.^{5,6,9,11–13}

Aβ_16–22 at N_pep = 16 was self-assembled in the presence of a series of short heparin chains (dp2, dp4, and dp8), to test if limiting heparin–peptide interactions would lead to commensurate effects on Q and beta-domain mixing, compared to dp18 heparin. Compared to baseline values for Aβ_16–22 in water, N_beta varied by less than ±1 (Fig. S12a, ESI†), but populations of beta-domains with Q ≥ 0.75 rose in proportion with heparin dp (Fig. 3(a)). Simultaneously, H_beta and beta-domain lifetimes rose in proportion with heparin dp (orange marks in Fig. 3(c)), thus confirming our proposed mechanism that heparin-peptide interactions prolong beta-domain lifetimes, allowing them to sample highly ordered structures.

Fig. 3 Distributions of Q for Aβ_16–22 at N_pep = 16 at varying heparin dp (a) and with rigid heparin analogues (b), represented by kernel density estimates (KDE). Scatter plots of mean H_beta against mean lifetimes of beta-domains for Aβ_16–22 at N_pep = 16 at varying heparin lengths (orange • annotated by dp) and with rigid analogues of dp18 heparin (green ■ annotated by *X, where X is the rigidity factor). The pink region in (c) indicates mean domain lifetime ±2 SE for N_pep = 16Aβ_16–22 in water.

N _comb also decreased as heparin's dp, H_beta, populations of domains with Q ≥ 0.75, and domain lifetimes increased, confirming the thermodynamic interpretation that heparin suppresses the S^mix of the beta-domains (Fig. S13a, ESI†). Heparin dp-dependent increases in Q and lifetime potentially signal proportionate fibril-enhancement, which has been reported in previous experimental studies of amylin and PACAP27 peptide fibrillation with heparins of varying lengths.^9,13

Next, we aimed to understand how heparin's rigidity influenced its effects on ordered Aβ_16–22 aggregation. This analysis was motivated by the knowledge that chitosan, a charged polysaccharide of greater rigidity than heparin (persistence lengths of chitosan and heparin are 6 nm⁴² and 4.5 nm,⁴³ respectively), is a strong inhibitor of Aβ fibrillation^5,16—a contrast to heparin's fibril enhancing properties. In particular, we hypothesized that increasing heparin's rigidity would suppress its ability to spawn beta-domains with Q ≥ 0.75 in Aβ_16–22 oligomers.

To test this hypothesis, we self-assembled Aβ_16–22 at N_pep = 16 with a series of rigid dp18 heparin analogues—heparin₅₀, heparin₁₀₀, heparin₂₀₀ and heparin₃₀₀—created by increasing the force constants of the glycosidic backbone angles by “rigidity factors” of 50, 100, 200, and 300, respectively. While these rigid heparins were increasingly biased towards extended conformations, their morphologies were not affected in strict proportion to their rigidity factors. Rather, the most prominent effect on the rigid heparins was the exclusion of collapsed conformations that were common in the reference heparin (henceforth referred to as heparin₁), where heparin's radius of gyration, R_gyr, was between 1.6–1.8 nm (Fig. 4(a)).

Fig. 4 (a) Kernel density distributions of R_gyr of heparins at varying rigidity factors, over four independent 3 μs long trajectories. (b) Compactness of Aβ_16–22 oligomers at N_pep = 16 with different heparin rigidities (orange •). Pink regions denote the compactness of peptide aggregates in water. Insets show snapshots of oligomers (pink volumes) with compactness 0.8 and 0.7 with heparin₁ and heparin₂₀₀ (orange strands), respectively.

Although differences in N_beta were within ±1 beta-strand (Fig. S12b, ESI†), populations of Q ≥ 0.75 beta-domains with the rigid heparins were indeed lower than that with dp18 heparin₁ (Fig. 3(b)). Moreover, the populations of beta-domains in the 0.4 ≤ Q ≤ 0.75 range were much higher in the presence of the rigid heparins, particularly heparin₃₀₀, compared to Aβ_16–22 in water. Thus, increasing heparin's rigidity not only suppresses its ability to spawn highly-ordered beta-domains (Q ≥ 0.75) but also increases the likelihood of poorly ordered domains (0.4 ≤ Q ≤ 0.75) in Aβ_16–22 oligomers.

As for lifetimes, beta-domains formed with rigid heparins were intermediate between those formed in water and in the presence of heparin₁ (green marks in Fig. 3(c)). Similarly, beta-domains’ S^mix with the rigid heparins was lower than with heparin₁, and higher than in water (Fig. S13b, ESI†). Altogether, rigid heparins spawned oligomers containing poorly ordered and short-lived beta-domains, which could serve as poor seeds that inhibit fibrillation at long timescales, as seen with chitosan.^5,16 These results also conform to our previous work where chitosan was shown to suppress beta-strand counts and extended conformations—both of which count towards the order parameter Q—among Aβ_16–22 peptides.¹⁷

Interestingly, beta-domains with rigid heparins had about as many heparin contacts, indicated by overlapping error bars in Fig. 3(c), and similar distributions of contacts across bead types (Fig. S8 and S9, ESI†) as those with heparin₁. Yet, beta-domain lifetimes were lower, as if the rigid heparins were dp8 in length (Fig. 3(c)). To explain this discrepancy we looked at the shapes of the Aβ_16–22 at different heparin rigidities and considered their implications.

In Fig. 4(b), shapes of Aβ_16–22 oligomers are described by their compactness, i.e., the ratios of the smallest and largest moments of inertia where 1 corresponds to a spherical form and 0 corresponds to a rod-like form. Snapshots of oligomers in Fig. 4(b) serve as visual references for the shapes of oligomers at compactness values 0.8 and 0.7. When heparin chains were collapsible (heparin₁ in Fig. 4(a)), peptide–heparin contacts were maximized (Fig. 3(b)) while slightly promoting the spherical character of the Aβ_16–22 oligomers relative to Aβ_16–22 in water. This enhancement in compactness implies a corresponding reduction in hydrophobic solvent accessible surface area. On the other hand, in trying to maximize contacts with extended and non-collapsible heparin chains (heparin₅₀, heparin₁₀₀, heparin₂₀₀ and heparin₃₀₀ in Fig. 4(a)), Aβ_16–22 oligomers also adopted more extended rod-like shapes. The promotion of rod-like oligomers implied a corresponding increase in hydrophobic solvent accessible surface area. Thus, we reason that competition between heparin–peptide interactions and the hydrophobic effect would lead to frustration in the Aβ_16–22 oligomer, which could be responsible for the lower beta-domain lifetimes and associated shifts in the ensemble from highly ordered (Q ≥ 0.75) to poorly ordered (0.40 ≤ Q ≤ 0.75) structures.

Finally, we characterised how heparin's rigidity affected the early kinetic pathways of Aβ_16–22 self-assembly leading up to the consolidation of peptides into a single oligomer. The motivation for this analysis was the jump in the time taken for consolidation, τ_agg, from tens of nanoseconds in heparin₁ to hundreds of nanoseconds among its rigid analogues (Table 1).

Table 1 Characteristic aggregation time at increasing heparin rigidities, τ_agg in nanoseconds

Rigidity factor	τ agg (ns)
1	50 ± 20
50	200 ± 100
100	400 ± 200
200	400 ± 200
300	600 ± 400

---

The consolidation of peptides into a single aggregate is achieved via a combination of two modes: the heparin-independent mode, where peptides consolidate in the bulk without heparin's involvement, and the heparin-dependent mode, where peptides consolidate on heparin's surface. To illustrate the role of rigidity in the heparin-dependent mode of consolidation, we focus on two simulations of Aβ_16–22 aggregation at N_pep = 16: one with heparin₁ (Fig. 5(a), (c) and (e)), and another with heparin₂₀₀ (Fig. 5(b), (d) and (f)). In both cases, unconsolidated peptide aggregates condensed at distant sites along the length of the heparin chain within the first 25 ns.

Fig. 5 (a) and (b) Snapshots of Aβ_16–22 aggregates (red, purple or blue) bound to heparin (orange chain) at representative time-points. (c) and (d) Heparin subunit indices bound by peptide aggregates, marked in their respective colours from their snapshots in (a) and (b), over time. (e), (f) R_gyr of heparin molecules over time, concurrent with (a) and (b). Data from a single trial of Aβ_16–22 aggregation at N_pep = 16 with dp18 heparin₁ is shown on the left (a), (c) and (e). Data from a single trial of Aβ_16–22 aggregation at N_pep = 16 with dp18 heparin₂₀₀ is shown on the right (b), (d) and (f).

Each peptide aggregate, identified by its constituent peptide id's, was assigned a unique colour, as in the snapshots in Fig. 5(a) and (b) depicting key moments in the consolidation process. Concurrently, the indices of the heparin subunits in contact with each peptide aggregate were marked with the aggregate's assigned colour, as a time series (Fig. 5(c) and (d)). A third plot concurrently tracked the associated heparin's R_gyr as a function of time (Fig. 5(d) and (e)).

In heparin₁, there were three aggregates at time = 10 ns (brown, olive, and grey), which merged into two (red and purple) and finally consolidated into a single aggregate (blue) around time = 70 ns (Fig. 5(a) and (c)). At the same time, heparin₁ underwent a gradual collapse from time = 0 ns to time = 150 ns (Fig. 5(e)). Thus, the bending of heparin seemed to be the dominant force for peptide consolidation. In contrast, heparin₂₀₀ maintained an extended conformation throughout (Fig. 5(f)) ruling out any contribution from its bending motions. Instead, a consolidated aggregate (blue) formed around the 400 ns mark exclusively by the larger purple aggregate crawling along the heparin₂₀₀ chain towards the red aggregate (Fig. 5(b) and (d)). These aggregation events with heparin₁ and heparin₂₀₀ were captured in Movies S1 and S2 (ESI†), respectively.

In effect, the heparin-dependent consolidation is composed of two pathways: one where heparin bends to merge aggregates at distant sites along the heparin chain, and another where aggregates crawl along heparin chains to merge. Heparin index occupancy and R_gyr data for all trials with all heparins (ESI,† Fig. S16–S20) demonstrate that there is a mix of both pathways where heparin is flexible (heparin₁), but the balance between the two pathways shifts towards peptide crawling at higher rigidity factors. Even in Fig. 5(c), there is evidence of the red and purple aggregates crawling along heparin₁ between 30 and 60 ns, albeit somewhat obscured by heparin₁ 's dramatic collapse.

Thus, the delay in the consolidation of peptides into a single aggregate among rigid heparins (Table 1) was attributed to the inhibition of their bending motions and the associated reliance on aggregate-crawling as the pathway for consolidation. Looking back, we can associate the early kinetic pathways among rigid heparins with less ordered (Fig. 3) and more frustrated and rod-like oligomers (Fig. 4(b)). In other words, heparin's ability to enhance ordered Aβ_16–22 aggregation is severely compromised without its flexibility.

4 Conclusions

In conclusion, we describe a potential mechanism underlying heparin's fibril-enhancing activities on Aβ,^5,6 and other amyloidogenic peptides.^12,13,44 We showed that contacts of Aβ_16–22 peptides with heparin promoted long-lived and highly ordered beta-domains, i.e., locally ordered beta-sheet-like regions within oligomers. A thermodynamic interpretation was also outlined, whereby heparin suppressed the entropy of mixing, S^mix of the beta-domains within the peptide oligomers. By promoting ordered oligomers, heparin could spawn better seeds for fibrillation relative to peptide aggregation in water.

Additionally, results from our tests of Aβ_16–22 aggregation with rigid heparin analogues may help us understand why chitosan is an inhibitor of Aβ aggregation.^5,16 As with the rigid heparins, rigid chitosan chains would increase the peptides’ reliance on crawling, thus delaying their consolidation, ultimately resulting in frustrated and less ordered aggregates. Furthermore, chitosan tends to self-assemble into hydrogel networks. The complex topologies of hydrogel networks could exacerbate the delays in aggregation to such an extent that peptide aggregates could become quasi-sequestered, which we observed in our previous simulations of Aβ_16–22 aggregation with chitosan.¹⁷

As mentioned in the introduction, the two opposite effects, of fibril enhancement by heparin and fibril inhibition by chitosan, have been demonstrated in several proteins including Aβ₄₂,^5,16 tau_298–317,⁷ α-synuclein⁸ and amylin,¹⁸ whose net-charges range from −9 for α-synuclein to +6 for amylin. In light of these peptides’ diversity in sequence and net-charge, we posit that the mechanisms of polysaccharide rigidity-dependent peptide aggregation demonstrated here, with the net-charge 0 model peptide Aβ_16–22, may be independent of the protein sequence and net-charge.

The methods and insights generated in this paper may also help us understand and potentially harness heparin's more recently discovered abilities to induce specific fibril polymorphs and modulate liquid–liquid phase separation of proteins.^14,15

Author contributions

S. M. conceived and acquired funding for this research. S. G. performed investigations and drafted the manuscript. S. G. and S. M. co-edited the manuscript and performed data analyses. A. P. and I. B. performed investigations and data analyses towards heparin model development.

Data availability

Outputs of molecular dynamics simulations generated in this article are available on Zenodo at https://doi.org/10.5281/zenodo.11479853. Input files for the simulations and a brief tutorial are available at https://github.com/suhasgotla/heparin_amyloid_self-assembly.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This research was supported by NSF award CHE-2202281 and supercomputing resources provided by the University of Maryland (https://hpcc.umd.edu).

Notes and references

1. R. U. Margolis, R. W. Ledeen, M. Sbaschnig-Agler, M. C. Byrne, R. L. Klein, B. H. Douglas II and R. K. Margolis, J. Neurochem., 1987, 49, 1839–1844.
2. W. L. Kiang, T. Krusius, J. Finne, R. U. Margolis and R. K. Margolis, J. Biol. Chem., 1982, 257, 1651–1659.
3. J. Daz-Nido, F. Wandosell and J. Avila, Peptides, 2002, 23, 1323–1332.
4. H. Wang, P. Cao and D. P. Raleigh, J. Mol. Biol., 2013, 425, 492–505.
5. J. J. Valle-Delgado, M. Alfonso-Prieto, N. S. de Groot, S. Ventura, J. Samitier, C. Rovira and X. Fernàndez-Busquets, FASEB J., 2010, 24, 4250–4261.
6. B. Klajnert, M. Cortijo-Arellano, M. Bryszewska and J. Cladera, Biochem. Biophys. Res. Commun., 2006, 339, 577–582.
7. M. Islam, E. Argueta, E. P. Wojcikiewicz and D. Du, ACS Chem. Neurosci., 2022, 13, 3034–3043.
8. S. Mehra, D. Ghosh, R. Kumar, M. Mondal, L. G. Gadhe, S. Das, A. Anoop, N. N. Jha, R. S. Jacob, D. Chatterjee, S. Ray, N. Singh, A. Kumar and S. K. Maji, J. Biol. Chem., 2018, 293, 12975–12991.
9. S. Jha, S. M. Patil, J. Gibson, C. E. Nelson, N. N. Alder and A. T. Alexandrescu, J. Biol. Chem., 2011, 286, 22894–22904.
10. E. Bazar and R. Jelinek, ChemBioChem, 2010, 11, 1997–2002.
11. N. Quittot, M. Sebastiao and S. Bourgault, Biochem. Cell Biol., 2017, 95, 329–337.
12. S. K. Maji, M. H. Perrin, M. R. Sawaya, S. Jessberger, K. Vadodaria, R. A. Rissman, P. S. Singru, K. P. R. Nilsson, R. Simon, D. Schubert, D. Eisenberg, J. Rivier, P. Sawchenko, W. Vale and R. Riek, Science, 2009, 325, 328–332.
13. M. Sebastiao, N. Quittot, I. Marcotte and S. Bourgault, Biochemistry, 2019, 58, 1214–1225.
14. Y. Tao, Y. Sun, S. Lv, W. Xia, K. Zhao, Q. Xu, Q. Zhao, L. He, W. Le, Y. Wang, C. Liu and D. Li, Nat. Commun., 2022, 13, 4226.
15. D. K. Garg and R. Bhat, Biophys. J., 2022, 121, 2568–2582.
16. H. Liu, B. Ojha, C. Morris, M. Jiang, E. P. Wojcikiewicz, P. P. N. Rao and D. Du, Biomacromolecules, 2015, 16, 2363–2373.
17. S. Gotla and S. Matysiak, Phys. Chem. Chem. Phys., 2023, 25, 10113–10120.
18. Q.-Y. Meng, H. Wang, Z.-B. Cui, W.-G. Yu and X.-Z. Lu, Molecules, 2020, 25, 1314.
19. M. Schäffler, S. Samantray and B. Strodel, Int. J. Mol. Sci., 2023, 24, 11238.
20. S. Samantray and B. Strodel, J. Phys. Chem. B, 2021, 125, 5511–5525.
21. X. Dong, R. Qi, Q. Qiao, X. Li, F. Li, J. Wan, Q. Zhang and G. Wei, Phys. Chem. Chem. Phys., 2021, 23, 20406–20418.
22. B. Khurshid, A. U. Rehman, R. Luo, A. Khan, A. Wadood and J. Anwar, ACS Omega, 2022, 7, 15132–15144.
23. A. Sahoo, P.-Y. Lee and S. Matysiak, J. Chem. Theory Comput., 2022, 18, 5046–5055.
24. S. Khan, J. Gor, B. Mulloy and S. J. Perkins, J. Mol. Biol., 2010, 395, 504–521.
25. J. J. Balbach, Y. Ishii, O. N. Antzutkin, R. D. Leapman, N. W. Rizzo, F. Dyda, J. Reed and R. Tycko, Biochemistry, 2000, 39, 13748–13759.
26. X. Zhou, Y. Wang, W. Zheng, G. Deng, F. Wang and L. Jin, Front. Mol. Biosci, 2022, 9, 824146.
27. A. V. Nogueira, D. L. Drehmer, M. Iacomini, G. L. Sassaki and T. R. Cipriani, Carbohydr. Polym., 2017, 157, 72–78.
28. D. H. Atha, B. Coxon, V. Reipa and A. K. Gaigalas, J. Pharm. Sci., 1995, 84, 360–364.
29. S. O. Yesylevskyy, L. V. Schäfer, D. Sengupta and S. J. Marrink, PLoS Comput. Biol., 2010, 6, e1000810.
30. D. H. de Jong, G. Singh, W. F. D. Bennett, C. Arnarez, T. A. Wassenaar, L. V. Schäfer, X. Periole, D. P. Tieleman and S. J. Marrink, J. Chem. Theory Comput., 2013, 9, 687–697.
31. S. J. Ganesan and S. Matysiak, J. Chem. Theory Comput., 2014, 10, 2569–2576.
32. S. J. Ganesan, H. Xu and S. Matysiak, Phys. Chem. Chem. Phys., 2016, 18, 17836–17850.
33. S. Khan, K. W. Fung, E. Rodriguez, R. Patel, J. Gor, B. Mulloy and S. J. Perkins, J. Biol. Chem., 2013, 288, 27737–27751.
34. D. J. Evans and B. L. Holian, J. Chem. Phys., 1985, 83, 4069–4074.
35. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
36. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
37. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
38. S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess and E. Lindahl, Bioinformatics, 2013, 29, 845–854.
39. N. Michaud-Agrawal, E. J. Denning, T. B. Woolf and O. Beckstein, J. Comput. Chem., 2011, 32, 2319–2327.
40. P. H. Nguyen, M. S. Li, G. Stock, J. E. Straub and D. Thirumalai, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 111–116.
41. G. M. Castillo, W. Lukito, T. N. Wight and A. D. Snow, J. Neurochem., 1999, 72, 1681–1687.
42. G. Berth, H. Dautzenberg and M. G. Peter, Carbohydr. Polym., 1998, 36, 205–216.
43. G. Pavlov, S. Finet, K. Tatarenko, E. Korneeva and C. Ebel, Eur. Biophys. J., 2003, 32, 437–449.
44. N. N. Jha, A. Anoop, S. Ranganathan, G. M. Mohite, R. Padinhateeri and S. K. Maji, Biochemistry, 2013, 52, 8800–8810.

---

Footnote

† Electronic supplementary information (ESI) available: Heparin model development, supporting data on the properties of Aβ_16–22 aggregates and movies depicting aggregation pathways. See DOI: https://doi.org/10.1039/d4cp02331e

---