O-glycans Expand Lubricin and Attenuate Its Viscosity and Shear Thinning

Abstract

Lubricin, an intrinsically disordered glycoprotein, plays a pivotal role in facilitating smooth movement and ensuring the enduring functionality of synovial joints. The central domain of this protein serves as a source of this excellent lubrication and is characterized by its highly glycosylated, negatively charged, and disordered structure. However, the influence of O-glycans on the viscosity of lubricin remains unclear. In this study, we employ molecular dynamics simulations in the absence and presence of shear, along with continuum simulations, to elucidate the intricate interplay between O-glycans and lubricin and the impact of O-glycans on lubricin's conformational properties and viscosity. We found the presence of O-glycans to induce a more extended conformation in fragments of the disordered region of lubricin. These O-glycans contribute to a reduction in solution viscosity but at the same time weaken shear thinning at high shear rates, compared to nonglycosylated systems with the same density. This effect is attributed to the steric and electrostatic repulsion between the fragments, which prevents their conglomeration and structuring. Our computational study yields a mechanistic mechanism underlying previous experimental observations of lubricin and paves the way to a more rational understanding of its function in the synovial fluid.

Figure 1

Sequence, charge, and glycosylation of lubricin. (A) Schematic representation of the human synovial joint. Lubricin is found in the synovial fluid and acts as a lubricant to ensure the smooth movement of bones. (B) Lubricin is composed of two globular terminal domains (blue) and a mucin-like central disordered region (gray). The two globular domains are nonglycosylated. The central region spans approximately 800 amino acids, is highly glycosylated by negatively charged O-glycans, and is rich in polar and proline residues. Glycan sites, charged residues, and proline residues along the sequence are indicated. (C) Core I and core II of six O-glycans of lubricin (out of 11 types) used in this study. The NeuAC sugar has a negative charge (purple diamond). The net charge of each type of O-glycan and the prevalence for lubricin are indicated. (D) Five different 80-amino acid fragments spanning the disordered region of lubricin. They were glycosylated to a different extent (six levels of glycosylation for each fragment). The cumulative net charge for each glycosylated fragment is shown. (E) Examples of the initial fully extended structure for a nonglycosylated and a fully glycosylated fragment (sequence 461–540).

Figure 2

MD protocol to obtain viscosities of lubricin-derived fragments. (A) By using Monte Carlo sampling, we generated five different glycosylation distributions for each of the five 80-amino acid long lubricin fragments. Here, the blue line indicates the nonglycosylated fragment and the yellow symbols depict the O-glycans. (B) EMD simulations were carried out for each fragment, both in its nonglycosylated and its five different glycosylated forms. (C) EMD simulations were performed for systems containing multiple fragments in order to obtain the viscosity in the absence of shear by using the GK method. (D) Shear-driven NEMD simulation, deforming the simulation box, estimated viscosity under different shear rates.

Figure 3

O-glycosylation increases extension and stiffness of lubricin’s single-chain fragments. (A,B) Radius of gyration r_g as a function of the number of glycosylated residues N_g (A) and the net charge of the fragment q (B). Symbols display the values recovered from the MD simulations (average ± s.e., n = 3). Dashed lines correspond to a linear regression of the data of the form r_g = a_NN + b_N (A) and r_g = a_q|q| + b_q (B), with resulting fitting parameters a_N = 0.036 [nm] and b_N = 2.77 [nm] (A) and a_q = 0.036 [nm/e] and b_q = 2.72 [nm] (B). Cartoons exemplify compact and extended conformations for a fragment without and with bound sugars, respectively (protein: black, sugars: orange).

Figure 4

Viscosity calculation procedure at zero shear (A–C) and under shear flows (D,E). (A) Components (P_αβ) of the pressure tensor extracted from EMD simulations (shown here is the time-trace of one component for an exemplary system). (B) The average of autocorrelation of all pressure components is computed for each independent system. (C) The zero shear viscosity is obtained from the average autocorrelation of all pressure components using eq 3 (see Materials and Methods). Gray curves represent the independent viscosity of each replica (n = 50). The black curve displays the average of all these curves. The viscosity was extracted from the plateau-highlighted region (average ± standard error). Figure S7 shows all three cases of zero shear viscosity. (D) To estimate the shear viscosity, the simulation box was deformed at a shear rate . Initial (t = 0) and posterior (t > 0) snapshots are shown, highlighting the lubricin fragments in blue. (E) Velocity profiles were obtained upon box deformation under different shear rates. The line corresponds to the linear fit. Output shear rates, i.e., u/h, are shown in the legend. Also, the expected values of shear rates are 0.94, 0.625, 0.50, 0.25, 0.125, 0.06, and 0.03 ns^–1.

Figure 5

Shear thinning behavior of glycosylated lubricin fragments. (A) Viscosity η as a function of the shear rate obtained from equilibrium simulations and shear-driven NEMD simulations . Five different systems were considered, either without glycans (“w.o”: open symbols and dashed or dotted lines) or with glycans (“w.”: closed symbols and solid lines). Systems with medium (“w.”: green) and high (“w.+”: orange) levels of glycosylation were considered at the indicated mass (ρ) and molar densities (with N indicating the number of peptides in the system). A low-density system without glycans was also simulated (blue). For systems where zero shear Newtonian viscosities are available, fits to the Carreau model (eq 6, a = 2) are shown (solid or dotted lines). Dashed lines are fitted to a simple power law expression without a Newtonian plateau. The viscosity of water, both experimental and obtained here, is also displayed as reference (gray and black, respectively). (B) Sample geometry for Reynolds calculations with a flat wall sliding at velocity U against a parabolic height profile with h_min = 50 nm, h_max = 10 μm, and L_x = 10 mm. (C) Pressure profiles for lubricants with medium and high levels of glycosylation and three different sliding velocities. At low speed, normalized pressure profiles (by a reference pressure η₀LU/h²) fall onto the same curve, indicating that the flow is still in the Newtonian regime, as shown in the inset. (D) Effect of shear thinning in the Reynolds calculations shown by the local relative viscosity along the sliding direction.

Figure 6

Shear-dependent rheological properties of lubricin fragments derived from MD simulations. (A,B) Variation of the radius of gyration (R_g) (A) and the ratio between the SASA of whole proteins and the cumulative SASA of individual proteins (surface exposure ratio, see Materials and Methods section) (B) as a function of the shear rate . Low values of S_tot/∑S_i indicate a high degree of protein aggregation, whereas a value of one denotes zero aggregation. Symbols represent data obtained from EMD simulations at zero shear rate and shear-driven NEMD simulations at nonzero shear rates (average ± standard error, n = 4). Dashed lines correspond to nonglycosylated systems (“w.o.”), while solid lines represent glycosylated ones (medium glycosylated: “w.” and highly glycosylated: “w.+”). Color indicates mass densities. The viscosity–shear response (of Figure 5A) is shown at the top of A for comparison. (C,D) The ratio of viscosity η(w.)/η(w.o.) is presented as a function of the ratio X(w.)/X(w.o.), with X = R_g (C) and X = S_tot/∑S_i (D), in both the medium (green) and high (orange) mass density regimes. Correlation coefficients (r) for each data set are indicated. (E) Representative snapshots for the extreme cases highlighted with arrows in A and B are shown (backbone: blue and (un)glycosylated side chains: gray).

Figure 7

Shear-dependent alignment of lubricin fragments derived from MD simulations. (A) NCF for bulk systems of nonglycosylated (top row) and glycosylated (bottom row) systems with different mass concentrations (columns) as a function of the radial distance. Color indicates the different shear rates. An NCF value of 0.5 indicates random chain orientation, while 1.0 indicates full alignment. The gray area corresponds to nearest neighbors region for which high alignment was observed. Representative snapshots for the extreme cases, highest and lowest shear rates, are shown (backbone: blue and (un)glycosylated side chains: gray). (B,C) NCF for short-range (<4 nm) (B) and for long-range (>4 nm) (C) interchain separations as a function of the shear rate . Symbols depict data obtained from the simulations (average ± standard error, n = 4). Dashed lines correspond to nonglycosylated systems (“w.o.”), while solid lines represent glycosylated “w.” systems. Colors indicate mass density. The viscosity–shear response (of Figure 5A) is shown at the top of B for comparison. (D,E) The ratio η(w.)/η(w.o.) is presented as a function of the ratio X(w.)/X(w.o.), with X = NCF_max, i.e., the maximum of the NCF, at short-range (C) and long-range (D) separations, in both the medium (green) and high (orange) mass density regimes. Correlation coefficients (r) for each data set are indicated.