In macroscopic methods, the EG is generally treated as a continuum layer with a limited number of bio-features, and the simulation outcomes are the macroscopic properties of the endothelial dendritic structures. Computational fluid dynamics (CFD; Tu et al., 2019) and the finite element method (FEM; Zienkiewicz et al., 2005) are two commonly used techniques in macroscopic simulations, which originated from the fluid mechanics and solid mechanics communities, respectively. Generally, continuum-based conservation laws and macroscopic physical models for bio-features are solved to obtain flow profile over and/or through the EG layer (Wang, 2007; Lee et al., 2016, 2020), deformation/force of the anchoring cells or membranes (Tarbell and Shi, 2013; Dabagh et al., 2014), order of magnitude of force transmitted via the glycocalyx (Saez et al., 2019).

In numerical modelling, simplifications of the problems under concern are inevitable due to the complexity of EG structures. With a focus on different aspects of the EG features, mechanotransduction can be approached from different perspectives. For example, the endothelial surface is covered by a layer of interweaving, dynamic and brush-like GAG chains. Such a layer can be simplified as an array of side chains or a layer of porous structures. By assuming the EG layer as an organised array of elastic chains, with individual chains having one end fixed to the cell membrane, Secomb et al. (2001) were able to identify the role of EG in transmitting shear stress to the endothelial cell surface and conclude that the EG layer would mediate rapid flow fluctuations and attenuate shear stress. By treating the EG layer as a porous medium, Tarbell and Shi reported that the solid stress of the endothelial cell surface is one to two orders of magnitude higher than the surface fluid stress in the presence of EG layer, suggesting that the interstitial flow shear stress can be sensed by the EG layer (Tarbell and Shi, 2013). As such, the role of EG is to strengthen the mechanical signals. Combining these results, a possible pathway for force transmission via the EG layer is proposed, i.e., flow shear stress → EG → cell membrane, as discussed in Tarbell and Pahakis (2006).

Simplification based on experimental observation is imperative for achieving numerical results of acceptable accuracy in reasonable timeframes. Take geometric periodicity as an example: Researchers (Lee et al., 2020) have set up different periodic geometries for core proteins and GAGs and analysed the pressure gradient (an indicator used to quantify mechanotransduction in the research) differences among different periodic geometries. Results suggest that the periodicity type determines the pressure gradient within the EG layer: hexagonal periodicity of the core proteins leads to greater mechanotransduction, and rectangular periodicity of aggregated GAGs leads to the greatest conversion of pressure gradients to wall tractions. Indeed, a hexagonal periodicity is favoured by previous studies (Weinbaum et al., 2003; Zeng et al., 2013), and conclusion based on the hexagonal periodicity is sensible.

For macroscopic methods, a “complete” model with all the features of the system is definitely beneficial for accurate results from numerical simulations. On the other hand, the computational cost would increase as more details are involved in computational models. Such a complete model is intimately related with access to accurate structural information which in turn relies on the availability (and possibly invention) of cutting-edge microscopic observation equipment.

The lattice Boltzmann method (LBM) and the dissipative particle dynamics (DPD) are two widely used mesoscopic methodologies in mechnotransduction studies. Both methods are based on the kinetic theory of particle dynamics.

In the framework of LBM, the fluid is modelled as a set of fictitious particles, and such particles interact with each other on a regular discrete lattice mesh through a two-stage “collision and streaming” process. A systematic review of the LBM in flow studies can be found in Chen and Doolen (1998) and Li et al. (2016). In DPD models, the system of interest is simplified to a set of interacting beads, each bead representing a cluster of constitutive molecules or atoms. The governing equation in the DPD method is Newton’s Law of Motion, i.e., F = ma, where F is the force on the bead, m is the mass of the bead, and a is the acceleration of motion of the bead. The force term comprises three elements: a conservative force, a random force and a dissipative force (Liu et al., 2014). A detailed review of DPD can be found in Espanol and Warren (2017). The particle-based nature of both mesoscale methods facilitates the construction of complex structures at the fluid-structure interface. Thus, recent mesoscale simulation studies are mostly about interpreting the fluid-structure interactions between the blood flow and the complex endothelium surface structures (Yin and Zhang, 2013).

At the mesoscales, more features of the glycocalyx can be included, which enables detailed comparison with experiments. By coupling the LBM with immersed boundary methods, alongside structural solvers, researchers revealed a broad range of behaviours for slender structures, from a single flap in a periodic array, to a small finite array of flaps, and finally to a large finite array (Yin and Zhang, 2013; O’connor et al., 2016; O’connor and Revell, 2019). The findings suggest that the flow instability over the slender array depends on the natural frequencies of the flow and the surface structures (O’connor and Revell, 2019). As a simple extension, one may expect that the shedding or impairment of EG chains would change flow patterns, cause flow instability and further disturb mechanotransduction. The role of EG in regulating flow field was also demonstrated in a recent DPD study (Jiang et al., 2021).

Beyond the fluid-structure interactions between the blood flow and EG layer, mesoscopic methods have been applied to a wide spectrum of haemodynamic problems, like cell suspension and leukocyte adhesion over endothelium (Yan et al., 2019), the glycocalyx–endothelium–erythrocyte interaction in the microcirculation (Pontrelli et al., 2015; Jiang et al., 2021) and the formation of the cell-free layer (Deng et al., 2012).

Molecular dynamics (MD) is a typical microscopic simulation method that describes how atoms and molecules move and interact during a period of time. The governing equation of MD is also Newton’s Law of Motion as in the DPD method, but a wide range of interatomic and intermolecular forces can be included and often parameterised as force fields. For the theoretical and numerical details of the MD method, the readers are referred to Frenkel and Smit (2002) and Rapaport (2011). When the structural configurations of biomolecules are known and a set of force field parameters are properly chosen, the trajectories of the constitutive atoms can be calculated by MD and the behaviour of the biomolecular system is then fully characterised. In principle, all biophysical problems could be tackled by MD if the structural information at the atomic scale were known. The bottleneck preventing the broad application of MD is its high computational cost, which limits not only the number of observed individual molecules to be included but also the duration of observations (i.e., simulated physical time). Despite such constraints, MD has been increasingly employed in studying dynamics of EG and surrounding molecules, thanks to its unique ability to gain atomic level insight. In the classic vasoregulation signalling problem – in what manner does glycocalyx transmit force to the cytoskeleton from the flowing blood?, the glycocalyx was regarded as a rigid body and the force was transmitted by the local torque induced by the movement of the glycocalyx core protein. In recent MD simulations, Jiang et al. (2020b) elucidated the dynamics of Syndecan-4 (a typical type of transmembrane protein of the glycocalyx), and concluded that a scissor-like motion of the Syndecan-4 transmembrane dimer, together with its bending, is responsible for the force transmission and the force transmitted into the cytoplasm is of an order 10–100 pN. The MD model and force transmission mode is summarised in Figure 2. The same all-atom model for MD has also been applied to understand how the local lipid membrane around glycocalyx core proteins responds to the flow shear stress (Jiang et al., 2020a). By artificially switching off the charges of the core protein, the contributions of electrostatic and van der Waals interactions to cell membrane deformation induced by flow were quantified for the first time (Jiang et al., 2020a). The capability of isolating specific effects in MD demonstrates an important advantage of in silico experiments in exploring the mechanisms of mechanotransduction.

It is worthwhile to mention that fluid-structure interactions between blood flow and EG can also be interrogated by MD. When the current supercomputing platforms are employed, large-scale MD simulations are able to tackle flow problems in the proximity of the EG layer. As reported in recent studies, the blood flow changes the conformation of EG (Cruz-Chu et al., 2014; Pikoula et al., 2018), and the EG spatial configuration in turn modifies the blood flow profile (Jiang et al., 2017, 2018). This is in accordance with the clearly detrimental cycle between the endothelial dysfunction and the degradation of glycocalyx (Zhang et al., 2018).

Together with MD, the Monte Carlo (MC) method, based on the stochastic theory, is also able to capture microscopic features of bio-systems. A recent MC study has demonstrated that the glycocalyx promotes cooperative binding and clustering of adhesion receptors (Xu et al., 2016). The role of glycocalyx in nanocarrier-cell adhesion is also investigated by the MC method (Agrawal and Radhakrishnan, 2007).