Molecular mechanism of claudin-15 strand flexibility: A computational study

Tight junctions are belt-like, circumferential cell-to-cell contacts between epithelial cells that regulate the passage of water, ions, and small molecules through paracellular space (Tang and Goodenough, 2003; Anderson and Van Itallie, 2009; Furuse, 2010; Zihni et al, 2016; Odenwald and Turner, 2017). Tight junction proteins polymerize into a two-dimensional network of branching strands that connect two adjoining cells (Staehelin et al, 1969; Claude and Goodenough, 1973; Tsukita et al, 2001; Furuse, 2010; Krystofiak et al, 2019). Tight junction transmembrane proteins include tetra-spanning proteins such as claudins and occludin and other mono- and tri-transmembrane domain-containing proteins (Furuse et al, 1993, 1998a; Ando-Akatsuka et al, 1996; Martìn-Padura et al, 1998).

Claudins are one of the major components of tight junctions that play a key role in tight junction assembly and are responsible for selective paracellular permeability (Furuse et al, 1998b, 1999; Sasaki et al, 2003; Van Itallie and Anderson, 2006; Angelow et al, 2008; Anderson and Van Itallie, 2009; Furuse, 2010; Günzel and Yu, 2013; Tsukita et al, 2019; Piontek et al, 2020). Based on studies where claudins are expressed in fibroblasts, it has been shown that claudin strands are dynamic with structural flexibility to bend, arch, and form new branches on the scale of nanometers to micrometers (Sasaki et al., 2003; Van Itallie et al., 2017; Zhao et al., 2018). The structural flexibility of claudin strands is important for the maintenance of barrier function during changes in cellular morphology and tissue rearrangement (Marchiando et al., 2011; Higashi et al., 2016; Varadarajan et al., 2019). In fact, variations in strand morphology can influence tight junction barrier function (Claude and Goodenough, 1973; Claude, 1978; Furuse et al, 2001; Weber et al, 2010).

Claudins are transmembrane proteins with four transmembrane helices (TM1–TM4) that form a left-handed bundle and two extracellular segments (ECS1 and ECS2). The extracellular segments consist of an extracellular helix (ECH), five antiparallel β-strands, and a highly conserved consensus W-LW-C-C motif (Colegio et al., 2002; Yu et al., 2009; Suzuki et al., 2014; Günzel and Yu, 2013). Both ECS1 and ECS2 are essential for claudin assembly and charge selectivity of the pores (Li et al., 2013; Suzuki et al., 2014; Suzuki et al., 2015; Krause et al., 2015; Piontek et al., 2008; Piontek et al., 2017; Angelow and Yu, 2009; Anderson and Van Itallie, 2009; Rossa et al., 2014).

The crystallographic structure of mouse claudin-15, a channel-forming claudin in small intestines, shows the linear arrangement of claudins in a row (Suzuki et al., 2014). This linear arrangement was later suggested to represent one of the intermolecular contacts between claudins in tight junctions (Suzuki et al., 2015). The authors suggest that the ECS2 loop and the structurally disordered loop in ECS1 are vital for tight junction pore formation. However, the atomic structure of the ECS1 loop was not revealed. Based on cysteine crosslinking experiments, Suzuki et al. (2015) proposed a set of intermolecular cis-interactions mediated by the β-sheets of ECS1 (here termed “face-to-face”) that lead to a continuous arrangement of claudins into half-pipe structures (Suzuki et al., 2015). Side-by-side cis-interactions between claudins in the same membrane form the strand backbone. The paracellular pore structure is then completed through head-to-head trans-interactions between extracellular loops of claudins on adjoining cells. This model was verified in atomic simulations confirming head-to-head interactions between claudin loop domains in the adjoining cells and elaborating on the molecular basis for paracellular ion channel selectivity (Samanta et al., 2018; Alberini et al., 2017; Alberini et al., 2018; Fuladi et al., 2020). However, the proposed model by Suzuki et al. (2015) does not account for the different morphologies of tight junction strands observed under freeze-fracture electron microscopy, including strands with curved structures at the micrometer length scales, and cannot explain the protein–protein interactions necessary for strand flexibility. Longer strand models and dynamics studies are needed to identify the molecular basis for claudin-15 strand flexibility.

In this study, we extended our previous work on mouse claudin-15 (Samanta et al., 2018) to simulate claudin strands at scales previously unexplored. This model of the strands is based on the high-resolution structure of claudin-15 in a state not bound to a structurally perturbing enterotoxin (Suzuki et al., 2014) and our previous work on refining its paracellular ion channels in the membrane (Samanta et al., 2018). We carried out molecular dynamics (MD) simulations of claudin strands ranging from 30 to 225 nm in length, with 36 to 300 monomers, in two parallel lipid membranes. We investigated their mechanical properties at microsecond timescales. To achieve longer time scales, we simulated the strands at all-atom as well as hybrid resolution, where the protein is represented atomically and lipid and water molecules are coarse-grained (Han et al., 2010; Wan et al., 2011; Han and Schulten, 2012). The curvature of claudin-15 strands was successfully visualized in simulation trajectories. We calculated the persistence length of the strands and identified intermolecular interactions that lead to their flexibility. These findings suggest a novel mechanism for the bending flexibility of claudin strands that is consistent with experimentally observed curvatures of tight junction strands.

MD simulations were performed using the program NAMD (Phillips et al., 2005). All-atom simulations were carried out using CHARMM36 forcefield for proteins (MacKerell et al., 1998; MacKerell et al., 2004; Best et al., 2012), ions (Jorgensen et al., 1983), and phospholipids (Klauda et al., 2010) with the TIP3P water model (Jorgensen et al., 1983). These simulations were carried out with a time step of 1 fs and assuming periodic boundary conditions. Langevin dynamics with a friction coefficient of γ = 5 ps⁻¹ was used to keep the temperature constant at 300 K. The Langevin Nosé–Hoover method (Feller et al., 1995) was used to maintain the pressure at 1 atm in constant pressure simulations. Long-range electrostatic forces were calculated using the particle mesh Ewald method (Darden et al., 1993) with a grid spacing of at least 1 Å in each direction. The simulations used a time step of 1, 2, and 4 fs for bonded, short-range nonbonded, and long-range electrostatic interactions calculations, respectively. A 1–4 rule is applied to the calculation of nonbonded interactions. The neighboring list was updated every 20 steps. The electrostatic interactions were switched to zero between 10 and 12 Å, and the van der Waals interactions cutoff was set to 12 Å. The distance cutoff for searching the lists of interacting pairs was set to 13.5 Å. Additional restraints were applied by enforcing a harmonic potential with a force constant of 2.0 kcal.mol⁻¹.Å⁻², unless otherwise stated.

The hybrid resolution models were built and simulated using the PACE forcefield (Han et al., 2010; Wan et al., 2011; Han and Schulten, 2012). PACE forcefield models lipids and solvents consistent with the MARTINI forcefield (Marrink et al., 2007), while proteins are represented by a united-atom (UA) model, where heavy atoms and polar hydrogens are explicitly represented. The hybrid resolution models were simulated using a time step of 2 or 3 fs, with the neighbor list being updated every 10 steps, and the van der Waals interactions cutoff was set to 12 Å. The distance cutoff for searching the lists of interacting pairs is set to 14 Å. Periodic boundary conditions were applied in all three directions. Langevin dynamics with a friction coefficient of γ = 5 ps⁻¹ was used to keep the temperature constant at 300 K. The Langevin Nosé–Hoover method (Feller et al., 1995) was used to maintain the pressure at 1 atm in constant pressure simulations. The dielectric constant is set to 15 as in MARTINI simulations with nonpolarizable water (Marrink et al., 2007).

Five models of claudin-15 strands in two parallel lipid bilayers were constructed. The initial structures of claudin-15 strands were based on our refined model of three claudin-15 pores (12 claudin-15 monomers) in two parallel POPC lipid bilayers (Samanta et al., 2018). This model was made of the mouse claudin-15 structure with very high sequence similarity to human claudin-15 (Samanta et al., 2018). In this model, each membrane contains six claudin-15 monomers that are arranged linearly in an antiparallel double-row and form three half-pipe structures. Claudin pores are then formed by head-to-head trans-interactions between extracellular domains of claudins in two parallel lipid bilayers, similar to the pores shown in Fig. 1. The three-pore claudin-15 model, equilibrated in POPC bilayers (Samanta et al., 2018), was replicated to build longer claudin strands consisting of 8, 20, 44, and 74 pores (each pore consists of four monomers) that were 27, 63, 135, and 225 nm long, respectively. Four patches of hydrated lipid bilayers (50 Å × 50 Å) were added on the two sides of the claudin strand to separate the strand from its periodic image. The systems were all neutralized by 150 mM NaCl salt.

The all-atom models were then converted into a hybrid-resolution model based on the PACE forcefield (Han et al., 2010; Wan et al., 2011; Han and Schulten, 2012). The conversion is performed using the Python scripts provided by Han’s laboratory at (http://web.pkusz.edu.cn/han/) to convert protein, lipids, and ions into the hybrid model. The system was then solvated by coarse-grained MARTINI water molecules (Marrink et al., 2007) using VMD (Humphrey et al., 1996). The hybrid model systems consist of 200K-1.1 M particles.

The smallest system, i.e., a system consisting of 36 claudin-15 monomers (eight pores) was simulated in all-atom representation, as well. This system was equilibrated for 300 ns in a multistep process, as described below. In the first step, the protein backbone and lipid head groups were restrained harmonically with a force constant of 2 kcal.mol⁻¹.Å⁻², and the lipid tails were equilibrated for 200 ps at constant volume and temperature. In the next step, the lipid head groups were released and the system was equilibrated at constant pressure for 2 ns, after which, the two extracellular loops of claudins (residues 33–44 of ECS1 and residues 148–155 of ECS2) were released to move freely while the rest of the protein backbone was gradually released by decreasing the force constant to 1.0, 0.75, 0.5, and 0.25 kcal.mol⁻¹.Å⁻² over 45 ns. After releasing all restraints, the system was equilibrated freely for 250 ns at constant pressure.

Four hybrid model systems containing 8, 20, 44, and 74 pores (36, 84, 180, and 300 claudin-15 monomers) were equilibrated for 1.125 µs following a process similar to the equilibration of the all-atom system. Briefly, by harmonically restraining the protein backbone and lipid head groups with a force constant of 2 kcal.mol⁻¹.Å⁻², the lipid tails were equilibrated for 200 ps at constant volume and temperature. The lipid head groups were then released, and the systems were equilibrated for 2 ns at constant pressure and temperature. Finally, the two extracellular loops of claudins (residues 33–44 of ECS1 and residues 148–155 of ECS2) were released to move freely, while the rest of the protein backbone was gradually released by decreasing the force constant to 1.0, 0.75, 0.5, and 0.25 kcal.mol⁻¹.Å⁻² over 70 ns. These equilibration processes were performed at a time step of 2 fs. The simulations were then followed by 1.055 µs of simulation time at constant pressure using a time step of 2 fs for the first 140 ns and a time step of 3 fs for the last 915 ns.

By applying the worm-like chain approximation of linear polymers, the persistence length (l_p) of claudin-15 strands is directly calculated from thermal fluctuations of the strands in equilibrium MD simulation trajectories. To achieve this, the claudin strand is assumed to be a two-dimensional linear polymer in the plane of a membrane made of discrete segments. Each of the four pore-forming claudins are considered to be a segment located at position $r→i,$ where $r→i$ is the center of mass of the four claudins.

At each point on the strand $r→i,$ the tangent angle ϕi is determined as the angle corresponding to the tangent vector $τi→=(r→i+1−r→i−1)|r→i+1−r→i−1|.$ The persistence length is then calculated as the length of the strand over which the tangent angles become uncorrelated. Assuming a discrete strand made of segments with length δs, the contour length between any two points, i and j, on the strand is estimated to be s = |i − j|δs. We defined the contour fragment δs to be the distance between two segments $δsi=12(|r→i+1−r→i|+|r→i−r→i−1|)$ averaged over all segments during the simulation time $δs=〈δs¯i〉.$

The factor of 2 in the denominator is due to the two-dimensional nature of the strands. The persistence length is calculated for the five claudin-15 strand simulations as described above. The last 600 ns of the simulation trajectories are loaded at every 60 ps, and the center of mass for each of the four pore-forming claudins is recorded for each frame. The initial slopes of the decay of the logarithm of the cosine correlation functions are used to calculate the persistence length. The number of points included in the curve fitting is determined by an R-squared cut-off of 0.95.

Convergence of the persistence length calculations was assessed by calculating the persistence length of the four systems (30, 60, 135, and 225 nm) over 150-ns time blocks. The 900-ns trajectories were divided into six 150-nm short trajectories and the persistence length was calculated for each block. The average persistence length of the four systems dropped from ∼520 nm for the first 150 nm trajectory to 137 ± 42 nm over the last four time blocks (600 ns), indicating that microsecond-long trajectories are generally required for proper sampling of claudin strands.

The relative orientation angle of adjacent claudins in the lipid membrane (θ) is calculated with respect to the initial arrangement of claudins in the crystal structure (θ = 0). For each pair of neighboring claudins located at positions i and i + 1, a vector $v→(i,t)$ is defined as connecting the center of mass of the backbone of residue 84 on TM2 of the first claudin (i) to the center of mass of the backbone of residue 170 on TM4 of the second claudin (i + 1) at time t. For each claudin (i), the orientation angle θ is defined as the angle between $v→(i,t)$ and $v→(i,t=0),$ which corresponds to the crystallographic arrangement (Fig. 5). The last 600 ns of the simulation trajectory of the longest strand (225 nm with 300 claudins) was used for this analysis with a frequency of 60 ps, and the angles were reported for each of the 296 claudins after excluding the edges.

Convergence of the trajectories was assessed again by dividing the trajectories into 200-ns time blocks and calculating the distribution of the relative orientation angles over each 200-ns trajectory. The results indicate that the distributions remain relatively unchanged over the last 600 ns of the simulations, especially for the longer strands.

To characterize the side-by-side cis-interactions between adjacent claudins, a contact map of residues from neighboring claudins that are within 4.5 Å of each other is created at each frame and correlated with the relative orientation angles (θ). The last 15 ns of the simulation trajectory of the longest strand, i.e., the 225 nm strand (300 claudins) was used for this analysis. To be consistent with the other analyses, the trajectory was analyzed every 60 ps. Claudin pairs were clustered based on their relative orientation angle (θ) in three groups corresponding to −30° < θ < −10°, −10° < θ < +10°, and +10° < θ < +30°. Based on the residue–residue contact maps at each cluster, the probability of the dominant contacts was determined, and the most dominant contacts are presented. The contacts with <20% probability of occurrence are not reported. None of these contacts involved residues in the upper half of the claudins.

The average minimum distance between ECS2 and ECS1 loops of the monomers forming the head-to-head trans-interactions was calculated as the average of the minimum distance between backbone atoms of the opposing ECS2 loops (residues 148–154) or the minimum distance between backbone atoms of opposing ECS1 loops (residues 39–42) in 225-nm long strand trajectory. The last 90 ns of the trajectory was used for this calculation, and the average is reported over 144 claudin pairs, excluding the four pairs at the edges of the strand. The number of contacts between ECS2–ECS2 loops and ECS1–ECS1 loops was also calculated for backbone atoms of the corresponding residues that are within 4.0 Å of each other.

The correlation function between the shape of the strands in the double-membrane system was calculated over the last 90 ns of the 225 nm strand. At each point on the strand, the tangent vector to the strand in each membrane is defined as the vector connecting two points, $r→1$ and $r→2,$ on adjacent claudins in a row. The points $r→1$ and $r→2$ are defined as the center of transmembrane helices in each claudin monomer represented here as the center of mass of four central residues (residue 12 on TM1, residue 92 on TM2, residue 122 on TM3, and residue 178 on TM4). The tangent vector $v→=r→1−r→2$ is then calculated for each of the two neighboring claudins. The correlation between the tangent vectors in the two membranes is then estimated as the normalized dot product of the vectors $v→top.v→bottom/|v→top‖v→bottom|.$ A value of 1 means the vectors are perfectly aligned and a value of 0 indicates that the two vectors are orthogonal.

The lateral displacement of the two double-row strands in the two membranes was calculated as the RMS distance between two points connecting the ends of the ECS2 loop (residues 136 on TM3 and 158 on TM4) in two opposing claudins over 148 pairs over the last 90 ns of the 225-nm strand simulation trajectory.

We used MD simulations to study the stability and mechanical properties of claudin strands consisting of claudin-15 pores in two parallel lipid membranes. We constructed atomic models of claudin strands with lengths ranging from 30 to 225 nm consisting of 36 to 300 claudin-15 monomers. These models are based on our previously refined model of three claudin-15 pores with 12 monomers (Samanta et al., 2018), which was replicated to create longer claudin strands. The functional characteristics of claudin-15 pores in this model were previously verified, showing that it recapitulates experimentally defined size and charge-selectivity of claudin-15 channels (Samanta et al., 2018).

Our MD modeling successfully modeled ion transport through the claudin-15 pore, but it was not able to model the long-term interactions and dynamics among claudin-15 monomers. To address these questions, we simulated the systems using a hybrid-resolution model, PACE, that combines a united atom representation of proteins with a coarse-grained model of lipids and solvent (MARTINI) to achieve longer time scales and better sampling of the configuration space, which is necessary to study strand organization in the membrane (Han et al., 2010; Wan et al., 2011; Han and Schulten, 2012). To validate the results of the PACE forcefield, the structural stability of claudin monomers, their fluctuations, and the integrity of claudin pores in the 36 claudin strands (30 nm long) are compared against a 250-ns all-atom trajectory of the same system. We then characterized the mechanical properties of the strands as a function of their length and determined intermolecular contacts between claudin monomers that confer flexibility to the strand.

A model of 36 claudins in two parallel POPC lipid membranes was constructed (Fig. 1 a) based on our previously refined three-pore model of claudin-15. In this model, each lipid bilayer contains an antiparallel double-row of claudin monomers arranged linearly through two sets of side-by-side cis- and face-to-face cis-interactions (Fig. 1 b). Head-to-head (trans) interactions with claudins in the second lipid membrane then result in the formation of claudin channels (pores) in the space between the two membranes. The scaffold of each pore in this model is made of four claudins, two from each membrane. The four monomers form the β-barrel-like structure of the channel. The ECS1 loops of four other monomers protrude into the channel and form head-to-head contacts between claudins of two membranes to help seal the ion conduction pore through hydrophobic interactions.

The system of 36 claudins was equilibrated for 250 ns with an all-atom representation of the system including proteins, lipids, water, and ions (Fig. 1 c). The system remained stable with an average root-mean-square deviation of 2.19 ± 0.49 Å with respect to the initial structure of the backbone (Fig. 2 a). Face-to-face cis-interactions between claudins in each membrane were maintained through hydrogen bonds between β-sheets (ECS1) of neighboring claudins with an average of 1.3 ± 0.5 hydrogen bonds. There was a least one hydrogen bond between the monomers ∼78% of the time. The stability of the system was comparable with the smaller systems of claudin strands (Samanta et al., 2018; Alberini et al., 2017).

In addition to this all-atom simulation, the 36-claudin system was simulated in a hybrid environment where the lipids and solvent molecules were coarse-grained and the protein was represented atomically (Han et al., 2010; Wan et al., 2011; Han and Schulten, 2012; Fig. 1 d). During this 1 µs-long simulation, the claudin strand remained stable with an average root-mean-square deviation of 4.21 ± 0.88 Å for the protein backbone with respect to the initial structure (Fig. 2 b). Further comparison of the two simulation trajectories indicates that the secondary structure of claudin monomers is well-preserved in the hybrid-resolution model (Fig. 2 c). The tetrameric structure of claudin pores is also maintained in both simulations, with average backbone root-mean-square deviations of 3.86 ± 0.30 and 6.56 ± 0.89 Å for all-atom and hybrid resolution trajectories after fitting the tetrameric claudin channels to the initial structure (Fig. 2, d and e).

Claudins form on channels that run parallel to the membrane surface and seal the paracellular space between two membranes to water and ions. A snapshot of the tetrameric structure of claudin pores in our simulations is shown in Fig. 2 e. Simulations of the 36-claudin systems show fully hydrated channels formed by extracellular domains of claudins that are ∼40–50 Å long (Fig. 3 a). In claudin strands, a series of parallel ion channels line the paracellular space forming ion conduction pores that are fully hydrated. The average linear density of water molecules along the pores is 7.5 water molecules/Å corresponding to an effective radius of 7.2–9.2 Å. The difference between the effective radii of the pore in all-atom and hybrid-resolution models is due to the coarse-grained nature of water molecules in the latter, where four water molecules are represented as one molecule (Marrink et al., 2007).

Fig. 3 b shows the density profile of water molecules on a plane crossing through claudin pores in the middle of two membranes. This cross-sectional plane runs parallel to the two membranes and cuts through eight claudin pores. These eight water-filled pores run parallel to this plane and are separated by regions with no water presence, i.e., proteins (dark blue). The channels are filled with water molecules at a density close to the bulk and are separated by extracellular domains (ECS1 and ECS2) of claudin. Head-to-head (trans) interactions between ECS1 and ECS2 of claudins in two opposing membranes seal the space between the two membranes and limit the flow of water and ions to claudin pores. Hydrophobic interactions between variable regions of ECS1 (residues 39–42) or residues 149–150 of ECS2 are responsible for head-to-head interactions of claudins (Furuse et al, 1999; Lim et al, 2008b; Piontek et al, 2008, 2017; Piehl et al, 2010). These interactions are well-maintained during the simulations and form well-defined pathways for ions and water in paracellular space despite the dynamic nature of the strands. Dynamic curvature of the strands is more apparent in the hybrid-resolution trajectories, where claudin tetramers forming ion channels adapt to the local curvature of the strand and maintain the pore integrity through persistent head-to-head and face-to-face cis-interactions.

To characterize the mechanical properties of claudin strands, three additional systems of 84, 180, and 300 claudins corresponding to strands of 63, 135, and 225 nm in length were prepared. Each system consists of two double-rows of claudins in two parallel POPC lipid membranes forming 20, 44, or 74 ion channels (Fig. 4 a). Each system was simulated for 1 µs in the hybrid-resolution representation, and the stability of the system was assessed similarly to the smaller system. The persistence length and bending rigidity of claudin strands in the plane of the membrane were calculated from thermal fluctuations of the system at 300 K.

The persistence length of the strands l_p is an indication of the stiffness (or flexibility) of the strand in the plane of the membrane and is defined as the length over which the direction of the strand becomes uncorrelated. The persistence length of claudin strands is calculated from equilibrium simulation trajectories. Assuming that the correlation function between tangent vectors at two points on the strand decays exponentially as a function of their separation, the persistence length is calculated as the characteristic length scale of this exponential decay along the strand. The persistence length of claudin-15 strands was estimated from the equilibrium trajectories of four strands and corresponds to an average of l_p = 137 nm with a standard deviation of 47 nm among four systems (Fig. 4 b). This is consistent with a persistence length of 191 ± 184 nm estimated from freeze-fracture electron microscopy images of claudin-15 strands (Zhao et al., 2018).

It has to be noted that Zhao et al. (2018) estimated the persistence length (l_p) based on the end-to-end distance of strands in FFEM images of tight junctions. This approach has been used to estimate the l_p of actin filaments and is easily correlated with the observed curvature of the strands. It requires, however, that the length of the strands be comparable to or larger than l_p. We have calculated the persistence length directly from the thermal fluctuations of the strands. In this approach, which is used for more rigid filaments, the strands do not need to be longer than l_p. In fact, this approach is more suitable for semirigid filaments, such as microtubules that are generally shorter than their persistence length (Howard, 2001). In the case of claudin strands, our calculations show that the persistence length calculated from the 30-nm strand is comparable with the values obtained from 136 and 225 nm strands.

Based on the persistence length, l_p, of the claudin-15 strands obtained from equilibrium trajectories, we estimated the bending rigidity κ_b of the strands as 563 ± 193 pN.nm² (Fig. 4 b). To put this value into context, ∼0.19 k_BT energy is needed to bend a 10-µm long claudin-15 strand through an angle of 30° at 300 K. This bending rigidity is comparable with the rigidity of double-stranded DNA with a persistence length of 50 nm (Frontali et al., 1979; Borochov et al., 1981; Manning, 2006; Geggier et al., 2011; Brinkers et al., 2009; Smith and Bendich, 1990), which is four orders of magnitude smaller than the bending rigidity of force-transmitting microtubules (Gittes et al., 1993; Venier et al., 1994; Sept and MacKintosh, 2010; Kikumoto et al., 2006; Felgner et al., 1996) and one order of magnitude smaller than that of actin filaments with persistence lengths in the millimeter and micrometer range (Gittes et al., 1993; Chu and Voth, 2005; Ott et al., 1993). This is consistent with the dynamic nature of claudin strands as seen in live fluorescence imaging (Sasaki et al., 2003; Van Itallie et al., 2017; Zhao et al., 2018).

Coarse-grained (CG) and all-atom simulations are commonly used to calculate the mechanical properties of biological systems from MD trajectories. For example, the bending rigidity of POPC lipid bilayers is reported as 29–32 k_BT from all-atom simulations based on the CHARMM36 forcefield (Venable et al., 2015; Fowler et al., 2016) and 20–25 k_BT from coarse-grained simulations using the MARTINI forcefield (Bochicchio and Monticelli, 2016; Fowler et al., 2016). Both forcefields consistently reproduced experimentally reported values of the POPC bending rigidities which are reported to be in the range of 5.8–49 k_BT (Loftus et al., 2013; Jablin et al., 2014). The PACE forcefield is shown to reproduce the interactions between proteins, represented atomically, and coarse-grained lipid and water molecules very accurately (Han et al., 2010; Wan et al., 2011; Han and Schulten, 2012). Hence, it is expected that the hybrid resolution model used here reasonably produces the persistence length and bending rigidity of claudin-15 strands.

The bending flexibility of claudin strands in the membrane stems from the flexible arrangement of monomers within the membrane. Upon bending, adjacent claudins in each row deviate from their initial arrangement (Suzuki et al., 2015; Suzuki et al., 2014) and are displaced to adapt to the local curvature. To characterize the dynamic behavior of claudin strands, we looked at the relative orientation of claudin pairs as the local curvature of the strand changes. We calculated the relative orientation angle of claudins as they rotate with respect to each other to accommodate local changes (side-by-side cis-interactions; Fig. 5, a–c), as well as the distance between opposing claudin pairs in the double-row arrangement of claudins in each membrane (face-to-face cis-interactions; Fig. 5, d and e).

The distribution of relative orientation angle between claudin pairs in the 225 nm strand (300 claudins) shows a symmetric distribution of angles between −70° and +70° (Fig. 5 c). The peak at the center at 0° corresponds to the linear arrangement of claudins observed in the crystal structure of claudin-15 (Suzuki et al., 2014). At a temperature of 300 K, 95% of adjacent monomers adopt a relative orientation of −50°^◦ to +50° and 68% have a relative angle of −26° to +26°. The width of the orientation distribution indicates that local curvatures up to 0.29 nm⁻¹ corresponding to orientation angles of up to 50° are permissible in claudin-15 strands. These data suggest that the side-by-side cis-interactions, present in the initial model, are highly flexible and support the range of curvatures observed in this work.

We further characterized the atomic interactions between claudin pairs forming side-by-side cis-interactions from the simulation trajectories. Based on the relative orientation angle of claudins θ (Fig. 5, a–c), we grouped claudin pairs into three clusters: one cluster corresponding to the linear arrangement observed in the crystal structure (−10° < θ < +10°), another cluster corresponding to the dominant positive orientation angles (+10° < θ < +30°), and the last cluster corresponding to the dominant negative orientation angles (−30° < θ < −10°). The three clusters represent three dominant interfaces between neighboring claudin pairs. Contact maps between claudin pairs show the main residues involved in the side-by-side cis-interaction of claudins in each cluster (Fig. 6, a–c). At each of these clusters, three sets of interactions/regions are identified (red, purple, and orange in Fig. 6). These sets include interactions of the ECH; residues P66, S67, M68, L69, and L71, with three distinct regions: (1) ECS2/TM4; residues K155, E157, L158, and Y163 (ECH–ECS2; purple cluster), or (2) TM3; residues T143, F146, and F147 (ECH–TM3; orange cluster), or (3) ECS1; residues T33, V34, and H35 (ECH–ECS1; red cluster). The probability of observing each of these interactions varies over time as the orientation angle between claudins changes.

The most dominant interaction is between residues S67/M68 of ECH and E157 of ECS2 (region I) consistent with previous reports (Samanta et al., 2018; Zhao et al., 2018). E157 is a highly conserved residue in the claudin family (Lim et al., 2008a; Piontek et al., 2008; Piontek et al., 2017; Krause et al., 2015), and in this model, is located at the entrance of the pore (Samanta et al., 2018). Mutation of E157 is shown to disrupt strand formation (Zhao et al., 2018). A second dominant interaction is between M68 of ECH and F146/F147 of TM3 (region II), which is observed mostly at the zero or positive orientation angles (Fig. 6, b and c). These interactions are reported in the crystal structure of claudin-15 (Suzuki et al., 2014) as the main connection between neighboring claudins for maintaining the cis-interactions. Mutation of M68, F146, and F147 to smaller or charged residues is shown to disrupt strand formation (Suzuki et al., 2014). We also observed a third dominant interface between ECH and ECS1 in the simulation, which involves the variable region of ECS1 residues 33–35 (region III). This region was not resolved in the crystal structure of claudin-15 and was modeled in our claudin-15 pore (Samanta et al., 2018). This set of interactions is mostly observed at zero or negative orientation angles (Fig. 6, a–c). The β1–β2 loop on ECS1 (32–42) is hypothesized to form a strong hydrophobic seal in trans-interactions (Krause et al., 2008; Yu et al., 2009; Suzuki et al., 2014), and MD simulations have identified strong pair-wise hydrogen bonds between residues 39 and 42 (the other end of the loop) of opposing claudins (Samanta et al., 2018). Therefore, residues 33–35 are accessible for side-by-side cis-interaction with the adjacent claudins. Snapshots of the three interfaces and the residues involved are shown in Fig. 6, d–f.

These findings suggest that the ECH acts as a pivot point around which claudin monomers rearrange and rotate around a centrally positioned ECH–ECS2 interface with S67–E157 as the dominant interaction site. At positive curvatures, corresponding to positive orientation angles (Fig. 6 c), the frequency of contacts between ECH and ECS1 decreases, and more contacts are established between ECH and TM3. Similarly, as the local curvature of strands changes and the relative orientation angle between neighbors becomes negative, the interactions between ECH and ECS1 become more dominant (Fig. 6 a). The spatial arrangement of ECS1 and TM3 on the two sides of ECS2 allows the strands to bend while maintaining side-by-side cis-contacts. These three interfaces, providing a wide range of permissible orientation angles for side-by-side cis-interaction, are critical in conferring flexibility to claudin strands.

Residues on ECH are critical for strand assembly and function. The proposed mechanism here further elucidates the role of these residues in providing a pivot point that allows claudins to rearrange in the strand. This is independent of the secondary structure of ECH which in more recent structures of claudin-19 (Saitoh et al., 2015), claudin-4 (Vecchio et al., 2021), and claudin-9 (Vecchio and Stroud, 2019) was modeled as an unstructured loop.

In addition to the side-by-side cis-interactions between claudins, we characterized the face-to-face cis-interactions between claudins monomers in the membrane (Fig. 5 d). The double-row arrangement of claudins involves an antiparallel double-row of claudins in each membrane that is stabilized through a network of four to six hydrogen bonds between β-sheets of claudins (face-to-face cis-interactions; Samanta et al., 2018). These hydrogen bonds are responsible for maintaining the tetrameric pore structure and forming a well-defined permeation pathway. We calculated the average distance between the two β-strands of claudins over the simulation trajectories of the most dynamic strand with 300 monomers. The average distance between the two monomers over the last 90 ns of simulation was 2.51 ± 1.57 Å with a median of 1.86 Å, indicating that face-to-face cis-interactions are well-maintained (Fig. 5 e). Similar to the all-atom simulation of the 30-nm system (36 claudins), hydrogen bonds between claudin pairs were maintained for ∼67% of the time with an average of 1.2 hydrogen bonds per claudin pair. These data indicate that the two claudin rows in each membrane move together and explain why the local curvature of the two rows is correlated in the simulation trajectories.

These findings are recapitulated in a schematic representation (Fig. 7) showing the two-dimensional movement of claudin monomers in the membrane with respect to each other. The arrows indicate the directional movement of tetrameric claudin channels (only a pair is shown here for clarity) in the plane of the membrane. Claudin pairs, coupled to each other through face-to-face interactions, slide with respect to each other in the plane of the membrane to produce a local curvature in the strands. Flexible side-by-side cis-interactions between claudins are key in this sliding motion. Three sets of side-by-side interactions between ECH and ECS1, ECS2, or TM3 facilitate this movement. ECH, shown in green, acts as a pivot point in maintaining the side-by-side interactions, while claudins slide in the membrane with respect to each other. During this movement, the tetrameric structure of claudin ion channels is maintained by head-to-head trans-interaction of the claudins between two membranes.

Another key interaction in accommodating the sliding motion without breaking the head-to-head trans-interactions or affecting the pore integrity is the face-to-face cis-interaction. The β4–β4 interaction between the monomers is maintained by hydrogen bonds between the backbone atoms of residues 61–64. However, similar to mechanical proteins, the two β4 strands can slide with respect to each other and shift the hydrogen bonding partners without breaking the contact between the monomers. This is slightly reflected in the fluctuations in the number of hydrogen bonds (±1) in the face-to-face cis-interactions. During this sliding motion, the head-to-head interactions between extracellular loops (ECS1–ECS1 and ECS2–ECS2 loops) are well maintained.

Lateral fluctuations and the local curvature of claudin strands are highly correlated between the monomers that reside in the two membranes. In the longest strand simulated here (225 nm), the orientation of the strand defined by its tangent vector is highly correlated between the two membranes with an average correlation of ∼0.93. This indicates a strong correlation between the shape of the top and bottom layer double-row strands during the 1-µs trajectories. In addition, the lateral displacement of the double-row strands between the membranes quantified by the RMS distance between opposing monomers was 3.6 ± 0.8 Å, indicating the head-to-head trans-interactions do not dissociate during the simulation.

The correlated movement of the strand in the double membrane system is due to head-to-head trans-interaction of the claudins monomers with those in a second membrane, which preserves the tetrameric structure of claudin pores as well. These head-to-head interactions are highly maintained during the simulations. Two main interaction sites are the short ECS2 loop between TM3 and TM4 (residues 148–154) and the β1–β2 loop on ECS1 (residues 39–42) that interact with ECS2 and ECS1 loops of opposing claudins. The average minimum distance between ECS2–ECS2 loops and ECS1–ECS1 loops are 3.6 ± 0.8 and 3.1 ± 0.5 Å during the simulations, respectively. In addition to maintaining a minimum distance, there are on average 2.3 ± 1.2 Å contacts between the backbone atoms of opposing ECS2 loops and 2.8 ± 1.0 contacts between opposing ECS1 loops in the two membranes, indicating that the head-to-head contacts are well-maintained in the simulations.

The hybrid-resolution simulations significantly increase the time- and length-scales accessible in our simulations without sacrificing the atomic resolution of protein–protein interactions. To assess the validity of the proposed mechanism in all-atom simulations, we compared the thermal fluctuation of the 30-nm strand (36 claudins) in all-atom simulation with that of the 225 nm (300 claudins) in hybrid-resolution simulation. The relative orientation angle of neighboring claudins in the all-atom system is normally distributed around θ = 0, similar to the hybrid-resolution system. However, the width of the distribution is only 40° in the all-atom system as compared to 130° in the hybrid-resolution simulations. This is expected since the hybrid models are simulated four times longer (1 µs) compared to the all-atom system (250 ns).

We also determined the most frequent side-by-side cis-interactions in the all-atom system. The most dominant interactions are between residues ECH (residues 67–68) and ECS2 (residues 157–158), between ECH (residues 67–68) and TM3 (residues 146–147), and between ECH (residues 67–68) and ECS1 (residues 34–35). This is consistent with the middle panel of Fig. 6 b corresponding to orientation angles between −10° < θ < 10° in the hybrid-resolution model and consistent with lower fluctuations of the strand in the all-atom system.

We have used molecular dynamics simulations to investigate the mechanical flexibility of claudin-15 strands ranging from 30 to 225 nm in length. Claudin strands within tight junctions are dynamic and respond to external forces in epithelial cells by rearranging within the cell membrane and adapting their curvature. The paracellular barrier function is expected to remain intact during this process. We used our previously refined model of claudin-15 paracellular ion channels (Samanta et al., 2018) to simulate claudin strands with 36–300 monomers in two parallel lipid membranes. To present a less computational yet accurate model of claudin strand assemblies, we simulated the strands at all-atom as well as hybrid resolution, where the protein is represented atomically, and lipid and water molecules are coarse-grained. Our results indicate that strands of claudin-15 are dynamic in the plane of the membrane with a persistence length of ∼137 nm, consistent with experimental studies (Zhao et al., 2018).

Our results reveal that the functional flexibility of claudin-15 strands is due to the interplay of three sets of intermolecular interactions; the side-by-side cis-interactions between claudins in a row as observed in the crystal structure of claudin-15, the face-to-face cis-interactions between β-sheets of claudins in two adjacent rows within the same membrane, and the head-to-head trans-interactions between extracellular domains of claudins in two membranes. The side-by-side cis-interactions are the most flexible interaction network among the three. They provide rotational flexibility to neighboring claudins with a pivot point at the ECH. The face-to-face cis-and head-to-head trans-interactions are the least affected upon curving of the strands and are responsible for maintaining the tetrameric structure of the pore.

These findings suggest a novel mechanism for bending flexibility of claudins strands, in which claudins are assembled into interlocking tetrameric ion channels along the strand that slide with respect to each other as the strands curve over submicrometer length scales. The side-by-side cis-interactions provide the structural flexibility for this movement while limiting its range. Water and ion density profiles indicate that during this process, the ion channels remain intact and the paracellular seal is maintained. Thus, the barrier function is preserved.

These studies provide a mechanistic view on the functional flexibility of tight junction strands and the role of claudin–claudin interactions in maintaining epithelial barrier function. Further studies are needed to understand the large-scale rearrangement of tight junction networks accompanied by the formation of new tight junction branches, ideally with full atomistic detail. Most importantly, the role of more complex lipid bilayers in modulating strand flexibility should be addressed.

Crina M. Nimigean served as the editor.

The work of F. Khalili-Araghi, S. Fuladi, S. McGuinness, and C.R. Weber was supported by the National Science Foundation grant MCB-1846021. This research is part of the Frontera computing project at the Texas Advanced Computing Center. Frontera is made possible by the National Science Foundation award OAC- 1818253.

The authors declare no competing financial interests.

Author contributions: All authors contributed to the writing of the manuscript. F. Khalili-Araghi and S. Fuladi designed the research. S. Fuladi carried out the simulations and analyzed the data. S. McGuinness analyzed the data.