Molecular Dynamics Simulations of Claudin-10a and -10b Ion Channels: With Similar Architecture, Different Pore Linings Determine the Opposite Charge Selectivity

Claudin polymers constitute the tight junction (TJ) backbone that forms paracellular barriers, at least for bigger solutes. While some claudins also seal the barrier for small electrolytes, others form ion channels. For cation-selective claudin-15 and claudin-10b, structural models of channels embedded in homo-polymeric strands have been suggested. Here, we generated a model for the prototypic anion-selective claudin-10a channel. Based on previously established claudin-10b models, dodecamer homology models of claudin-10a embedded in two membranes were analyzed by molecular dynamics simulations. The results indicate that both claudin-10 isoforms share the same strand and channel architecture: Sidewise unsealed tetrameric pore scaffolds are interlocked with adjacent pores via the β1β2 loop of extracellular segment 1. This leads to TJ-like strands with claudin subunits arranged in four joined rows in two opposing membranes. Several but not all cis- and trans-interaction modes are indicated to be conserved among claudin-10a, -10b, and -15. However, pore-lining residues that differ between claudin-10a and -10b (i.e., R33/I35, A34/D36, K69/A71, N54/D56, H60/N62, R62/K64) result in opposite charge selectivity of channels. This was supported by electric field simulations for both claudins and is consistent with previous electrophysiological studies. In summary, for the first time, a structural and mechanistic model of complete and prototypic paracellular anion channels is provided. This improves understanding of epithelial paracellular transport.

Recently, we compared two variants of Suzuki's JDR architectural model for the two related cation channel-forming CLDN10b and -15: a tetrameric-locked pore barrel vs. octameric-interlocked pore barrels model [15].Molecular dynamics simulations of membraneembedded dodecamers indicated that CLDN10b and CLDN15 share the same architecture of TJ-strands: Octameric Interlocked pore Barrels (IB) in which sidewise-unsealed tetrameric pore scaffolds are interlocked with adjacent pores by the β1β2 loop of the extracellular segment (ECS) 1 (Figure S2B,C).On the one hand, this loop mediates hydrophobic clustering and, together with ECS2, cisand trans-interaction between subunits of adjacent tetrameric pore scaffolds, leading to the formation of polymeric strands.On the other hand, the β1β2 loop contributes to the lining of the pathway for ion conduction and, thus, to selectivity properties [15].
For paracellular anion channels, no similar oligomer modeling and MD simulations have been performed so far.Simulations and free energy calculations have been performed for CLDN4 [21,28,29].However, they were restricted to tetramers and did not consider embedment of channels into strands, and the contribution of CLDN4 to anion channels is debatable [2,30,31].In addition, a pathogenic CLDN5 mutant was reported to form an anion channel [32], and corresponding elegant MD simulations and free energy calculations supported this conclusion [33].However, the physiological function of CLDN5 is to form a barrier, in particular at the blood-brain barrier [1,7,32,34].Thus, we aimed here to analyze the prototypic and functionally well-characterized anion channel formed by CLDN10a [35][36][37].For this purpose, we took advantage of the high sequence similarity between CLDN10a and CLDN10b-analyzed previously.The two isoforms differ only in transmembrane helix 1 (TM1) and ECS1, and the differing charge selectivity is mediated by ECS1 [38], which still shows ~57% sequence similarity between both isoforms.
Here, we generated and tested a CLDN10a dodecamer (three-pore) model based on the previously generated CLDN10b model and the previously established workflow, including homology modeling, membrane embedment, distance constraints, and Desmond MD simulations.The results indicate that CLDN10a and -10b channels share a common interlocked pore barrels JDR architecture.Nevertheless, next to conserved interaction motifs, some inter-subunit interfaces differ partly.In addition, pore-lining residues differ between CLDN10a and -10b channels and result in opposite charge selectivity.

Generation of Dodecamer Model for CLDN10a
Previously, we generated the Iinterlocked pore Barrels (IB) model of CLDN10b channels based on experimental support [12,39], the CLDN15 JDR strand model [10], homology modeling, and MD simulations [15].This dodecamer model consists of three adjacent interlocked pore tetramers with the middle pore stabilized and completely closed sidewise by residues of eight subunits (chains) (Figure S2B,C).Here, we generated a CLDN10a homology model using the CLDN10b IB model as a template.The rationale was the assumption that CLDN10a and -10b share a similar overall architecture due to the high sequence similarity (80% in the template region, Figure S1A), a similar fold prediction (Figure S1B), and similar function (paracellular ion channel).
The CLDN10a protein dodecamer was embedded in two POPC lipid bilayers, and NaCl and water were then added.The system was relaxed and equilibrated with gradually releasing constraints similar to those established previously [15] and described in the Methods and Figure S3.Five variant MD simulation production runs (each for 100 ns) of the interlocked pore barrels (IB) model were performed and compared: IB-1 (Interlocked pore Barrels model 1): The assumed face-to-face interface was slightly supported by a force constant of 1 kcal mol −1 Å −2 on the backbone of cysteine 61 (C61) in β4-strand.All other atoms in the system could move freely.In addition, the initial IB-1 model was refined further to increase the interface-wise symmetry of the model.Different constraints during the production runs led to four variants of the refined model: IB-2-this model has the same C61 backbone constraints as IB-1; IB-2+lic-this model has additional constraints on the linear-cis interface (backbone and β-carbon of I67 and F68, 1 kcal mol −1 Å −2 ); IB-2+hc-this model has additional constraints on the backbone of helices (1 kcal mol −1 Å −2 ); and IB-2+ohc-this model has additional constraints on the backbone of helices of the four outer subunits, i.e., subunits that do not line the middle pore (1 kcal mol −1 Å −2 ).
The IB-2 model variant was selected as the reference model, and, as an example, the 100 ns snapshot of the corresponding MD simulation is shown in Figure 1E-G.For all simulated model variants (IB-1, IB-2, IB-2+lic, IB-2+hc, and IB-2+ohc), all twelve CLDN10a subunits were embedded well in the two opposing membranes (hydrophobic transmembrane helices were covered with acyl chains of intact (water-free) lipid bilayer throughout the simulations).The overall arrangement of the CLDN10a dodecamer was similar to the JDR-like arrangement of CLDN10b, including interlocked barrel-like pores and lipids trapped in the central region between the two claudins rows of one membrane (Figure 1A-G).Compared to the pre-equilibrated starting structures and to the CLDN10b model, all interface types (including face-to-face-cis, linear-cis, and ECS2-ECS2-trans interfaces) were largely maintained.Also similar to CLDN10b, in CLDN10a dodecamers, stable clusters of the hydrophobic tips of β1β2 loops from four subunits were formed between neighboring pore centers.The latter were lined by hydrophilic residues (Figure 1C,G).However, the pore-lining residues differed strongly between the two claudin-10 isoforms.Of note, individual interfaces varied in detail over time between the individual subunits and between the five simulation lines (IB-1, IB-2, IB-2+lic, IB-2+hc, and IB-2+ohc).Thus, we evaluated and compared the MD simulations of the CLDN10a model variants in detail.(F) turned view on the middle pore.Protein chains were embedded well in membranes, and the overall arrangement, including trapped lipids (orange) between the two claudin rows, was similar to that of CLDN10b.In the CLDN10a pore, mainly Cl − ions (green spheres) were present, whereas in that of CLDN10b, mainly Na + ions (magenta spheres) were present, fitting to opposite charge selectivity of the two different channels.(D) Schema, showing the contact of β1β2 loops (red and green) between trans-interacting subunits (other claudin parts gray), highlighting the kinked/flat orientation of the loops towards subunits in the same membrane.This results in interlocked loops of neighboring pores (interlocked barrels (IB) arrangement).CDLN10a and -10b models showed similar IB arrangements.(C,G) Close-up of dodecamer centers.On both sides of the middle pore, hydrophobic clusters were formed similarly for CLDN10a and -10b by V37/I38 (CLDN10a) and V39/I40 (CLDN10b) residues (shown as surfaces) of trans-and cis-interacting β1β2 loop tips from four subunits (blue, green, red, and black).In contrast, hydrophilic pore-lining residues (sticks) differed strongly between CLDN10a and -10b.See also Video S1.

RMSD Indicates Overall Stability of CLDN10a Dodecamer Models
The stability of the individual subunits and of the oligomer during the production run was evaluated by calculating the root-mean-square-deviation (RMSD) of the protein backbone.The RMSD for the entire dodecamer with respect to its initial structure increased steadily for IB-1, whereas it reached a plateau for all IB-2 model runs after ~90 ns (Figure S4A).It was lowest (~1.0 Å) for the strongest constrained one (IB-2+hc) and highest (~2.2 Å) for the run with constraints on the linear-cis interface (IB-2+lic).For IB-2 and IB-2+ohc, it showed an intermediate RMSD of ~1.7 Å.These relative constant values were in the range of CLDN10b dodecamer models (1.5-2.2Å, Figure S4A) and slightly lower than the RMSD for CLDN15 dodecamer simulations (<3.5 Å, [26]).This indicated the overall stability of the membrane-embedded CLDN10a oligomers, which was most relevant for the weakly constrained IB-2 model variant.
In addition, RMSD was calculated with respect to the initial structure of each dodecamer subunit (Figure S4B).The RMSD and, thus, structural change over time differed Protein chains were embedded well in membranes, and the overall arrangement, including trapped lipids (orange) between the two claudin rows, was similar to that of CLDN10b.In the CLDN10a pore, mainly Cl − ions (green spheres) were present, whereas in that of CLDN10b, mainly Na + ions (magenta spheres) were present, fitting to opposite charge selectivity of the two different channels.(D) Schema, showing the contact of β1β2 loops (red and green) between trans-interacting subunits (other claudin parts gray), highlighting the kinked/flat orientation of the loops towards subunits in the same membrane.This results in interlocked loops of neighboring pores (interlocked barrels (IB) arrangement).CDLN10a and -10b models showed similar IB arrangements.(C,G) Closeup of dodecamer centers.On both sides of the middle pore, hydrophobic clusters were formed similarly for CLDN10a and -10b by V37/I38 (CLDN10a) and V39/I40 (CLDN10b) residues (shown as surfaces) of transand cis-interacting β1β2 loop tips from four subunits (blue, green, red, and black).In contrast, hydrophilic pore-lining residues (sticks) differed strongly between CLDN10a and -10b.See also Video S1.

Evaluation of the CLDN10a Dodecamer Models 2.2.1. RMSD Indicates Overall Stability of CLDN10a Dodecamer Models
The stability of the individual subunits and of the oligomer during the production run was evaluated by calculating the root-mean-square-deviation (RMSD) of the protein backbone.The RMSD for the entire dodecamer with respect to its initial structure increased steadily for IB-1, whereas it reached a plateau for all IB-2 model runs after ~90 ns (Figure S4A).It was lowest (~1.0 Å) for the strongest constrained one (IB-2+hc) and highest (~2.2 Å) for the run with constraints on the linear-cis interface (IB-2+lic).For IB-2 and IB-2+ohc, it showed an intermediate RMSD of ~1.7 Å.These relative constant values were in the range of CLDN10b dodecamer models (1.5-2.2Å, Figure S4A) and slightly lower than the RMSD for CLDN15 dodecamer simulations (<3.5 Å, [26]).This indicated the overall stability of the membrane-embedded CLDN10a oligomers, which was most relevant for the weakly constrained IB-2 model variant.
In addition, RMSD was calculated with respect to the initial structure of each dodecamer subunit (Figure S4B).The RMSD and, thus, structural change over time differed between individual subunits, reflecting a certain heterogeneity within the dodecamers.For each model variant, the mean RMSD of all subunits reached a plateau after ~70 ns with values (≤2 Å) in a range similar to that obtained in CLDN15 polymer simulations [27].The RMSD values of the transmembrane helices backbone with respect to the initial structure were slightly below the values of the whole protein backbone for each model variant, with mean RMSD values of ~1.5 Å throughout the last 30 ns for IB-1, IB-2, IB-2+lic, and IB-2+ohc and ~1.2 Å for IB-2+hc (Figure S4C).

Definition and Comparison of Interfaces within Dodecamer by Residue Contact Maps: cisand ECS2-trans-Interfaces
To define the interfaces between the subunits in the dodecamer, residue contact maps were generated for the regions of all subunits that are relevant to the JDR architecture model [10,15].
(i) First, the face-to-face-cis interface was analyzed.Along the pore pathway, the bottom and top of the pore are each formed by two claudin subunits connected by this antiparallel association of ECS1 β4-strands (Figures 2A and and S2A).At this interface, for CLDN10b, close distances (<5 Å) were found between the β4-strand residues S61 to D65 with the pair C63-C63 in the center (data taken from previous simulations [15], Figure 2B).Similarly, for CLDN10a, the corresponding residues F59 to P63 were close to each other, with the C61-C61 pair in the center.The mean distance for CLDN10a-IB-1 was slightly lower than for -IB-2 (Figure 2B).In sum, face-to-face cis-interfaces were similarly formed for the CLDN10a and -10b dodecamers.
(iv) Regarding β1β2 loop-ECS2 contacts, within one CLDN10a/-10b subunit, T32/T34, R33/I35 of β1β2 loop were close to E155/E157 in ECS2 for CLDN10a as well as for CLDN10b, indicating similar positioning of this β1β2 loop region relative to E155/E157 (Figure 2E).In trans, S35/G37, S36/T38 of β1β2 loop were close to E151/E153, Q152/Q154 in ECS2, similar to both CLDN10a and-10b.However, partly due to electrostatic attraction, The contact maps show a similar overall pattern for CLDN10b and -10a, reflecting a similar overall architecture of both channels.However, the pair distances are, on average, slightly higher for the CLDN10a dodecamer simulations IB-1 and IB-2.This and differences in contact details are evaluated below.
(i) For the face-to-face-cis interface type, the H-bonds between two β4-strands (F59 to P63 residue backbones) of subunit pairs of the middle pore were counted.The mean H-bond count per interface was ~ 1.7 for all models (Figure 3A) and thus in the same range as for CLDN10b and CLDN15 models (1-2, [15,27]).
(ii) For the linear-cis interface, the distance between F68 (Cγ, ECH region) and L156 (Cγ, ECS2 pocket) was measured over time.For all CLDN10a models, the distance was The contact maps show a similar overall pattern for CLDN10b and -10a, reflecting a similar overall architecture of both channels.However, the pair distances are, on average, slightly higher for the CLDN10a dodecamer simulations IB-1 and IB-2.This and differences in contact details are evaluated below.
(i) For the face-to-face-cis interface type, the H-bonds between two β4-strands (F59 to P63 residue backbones) of subunit pairs of the middle pore were counted.The mean

Definition and Comparison of Interfaces within Dodecamer by Residue Contact Maps: ECS1 Loop Clusters
The pore is in the central half on each side, bounded by cis-and trans-associations of β1β2-and β3β4 loops of four subunits (Figures 1G and 4A, [15]).Here, three different pairwise subunit combinations (i-iii) come into contact (Figure 4A,B-D headlines).For each combination, respective residue contact maps were generated for the CLDN10a IB-1 and IB-2 simulations and compared with those of previous CLDN10b simulations [15]: (i) Linear-trans pairs in a claudin row (Figure 4B): For CLDN10a/-10b residues in β1β2 loop, contact (<5 Å) of V37/V39, I38/I40 with V37/V39, I38/I40 and S36/T38, V37/V39 with W42/T44 and of V37/V39, I38/I40 with G53/T55, N54/D56 in β3β4 loop was detected.The pairwise distances were, on average, slightly higher for CLDN10a than for -10b.(ii) For the linear-cis interface, the distance between F68 (Cγ, ECH region) and L156 (Cγ, ECS2 pocket) was measured over time.For all CLDN10a models, the distance was constant over time, and for IB-2, IB-2+lic, IB-2+hc, IB-2+ohc~7.7 Å, and IB-1~6.6Å (mean of all interfaces in dodecamer, Figure 3B).However, the mean distance differed between individual linear-cis interfaces within a dodecamer, as shown for IB-2 as an example (Figure 3C).In addition, individual interfaces with at least one electrostatic interaction between the ECS2 pocket (E155) and ECH region (T66, I67, F68) were counted over time.For the different models, 4 ±1 out of 8 possible counts were obtained mostly and relatively constant over time (Figure 3D).Together, the measurements indicated that the linear-cis interface was largely maintained during the simulations for all CLDN10a model variants, similar to those shown for the corresponding CLDN10b model previously [15].
(iv) To evaluate the combined hydrophobic interactions mediated by the ECS2 (linearcis with ECH region and trans with opposing ECS2), we calculated the solvent accessible surface area (SASA) of interfacial residues (measure for water exclusion at the interface, Figure 3F).SASA (normalized to the side chain atom numbers) was lowest for the hydrophobic pocket (F144, 156).Low SASA was also obtained for ECH region residues I67 and F68, whereas V70 had much higher SASA.Trans-interfacial P147, F149, and V150 showed intermediated SASA with the lowest values for P147.As comparisons, (a) the corresponding non-interacting residues ("free" in Figure 3F) in the peripheral subunits and (b) the close-by positively charged K69 (ECH region) and R33 (β1β2 loop) had much higher SASA.Comparing the different CLDN10a model variants, the more constrained models (IB-2+lic, IB-2-hc, and IB-2-ohc) had lower sum SASA than the weakly constrained models IB-1 and IB-2 (Figure 3F).The data indicate that I67, F68, F144, F145, P147, and L156 strongly contribute to hydrophobic cis-/transinteraction in the ECS2 region.

Further Analysis of β1β2 Loop Clusters
For CLDN10b, it was suggested that the β1β2 loop tips (conserved among classic claudins) of four subunits form a hydrophobic cluster [15].Thus, the contribution of the β1β2 loop tip (V37, I38) of CLDN10a to hydrophobic cluster formation was analyzed.As mentioned above, the four loop tips can come into contact in three different combinations (Figure 4).As an example, for the crosswise-cis contact, the mean pairwise distances were calculated for V37-V37, I38-I38, V37-I38, I38-V37 (Cβ-atoms) and the different CLDN10a model variants (Figure 5A-C).The distances were between 7.2 and 14.3 Å.For all CLDN10a models, the I38-I38 distance was lowest (7.2-8.2Å), and V37-V37 distances were highest (11.4-14.3Å).The pairwise distances were similar to that of CLDN10b [15].
As a measure for hydrophobic clustering of the four loop tips, water exclusion was analyzed by SASA measurements.The mean SASA (during simulation time) for V37 and I38 of four interacting subunits were calculated for the two clusters in the dodecamer (left and right side of the middle pore, Figure 1G).For all models, the mean SASA differed considerably between the two clusters (42-74 Å 2 and 110-127 Å 2 , Figure 5D).The lower value was similar to the values obtained for the corresponding CLDN10b models [15].The higher value was still much lower than the SASA obtained for alternative CLDN10b models with a different β1β2 loop arrangement (195-411 Å 2 ), [15]).However, the differing SASA for the two clusters of one dodecamer indicates an unexpected asymmetry.The latter might be due to the limited accuracy of the simulations.Thus, CLDN10a IB-2, which showed the smallest difference between the SASA of the two clusters, was chosen as the reference model.
In sum, the results indicate the formation of similar hydrophobic β1β2 loop clusters for CLDN10a and CLDN10b, supporting the contribution of these conserved interactions to oligomerization.As a measure for hydrophobic clustering of the four loop tips, water exclusion was analyzed by SASA measurements.The mean SASA (during simulation time) for V37 and I38 of four interacting subunits were calculated for the two clusters in the dodecamer (left and right side of the middle pore, Figure 1G).For all models, the mean SASA differed considerably between the two clusters (42-74 Å 2 and 110-127 Å 2 , Figure 5D).The lower value was similar to the values obtained for the corresponding CLDN10b models [15].The higher value was still much lower than the SASA obtained for alternative CLDN10b models with a different β1β2 loop arrangement (195-411 Å 2 ), [15]).However, the differing SASA for the two clusters of one dodecamer indicates an unexpected asymmetry.The latter might be due to the limited accuracy of the simulations.Thus, CLDN10a IB-2, which showed the smallest difference between the SASA of the two clusters, was chosen as the reference model.
In sum, the results indicate the formation of similar hydrophobic β1β2 loop clusters for CLDN10a and CLDN10b, supporting the contribution of these conserved interactions to oligomerization.

Ion Permeation Pathway of Pore in CLDN10a Dodecamer Models
The ion permeation path in the CLDN10a dodecamer models (IB-2, IB-2+lic, IB-2+hc, and IB-2+ohc) was inspected by analyzing the diameter along the axis of the middle pore in the last 50 ns of the simulation (using HOLE).The pore profiles of the models differed up to 2.5 Å with respect to the mean diameter at particular positions along the pore axis (Figure 6B).The mean diameter differed also on both sides of the pore center that was surrounded by four N54 and four H60 residues (Figure 6B-D).For the reference model IB-2, the center showed a mean diameter of 6.7 Å.For IB-2 and IB-2+hc, on one side, the R62 region formed the narrowest site (5.1 Å for IB-2), whereas on the other side, the R62 region was wider (7.4 Å for IB-2) than around H60. Nevertheless, for all models, the mean pore diameter in the inner pore region (spanning R33-R62-H60-R62-R33) was smaller than in the pore periphery beyond K69 residues.To what extent do the mentioned variations R33 residues belong not to the tetrameric pore scaffold but to neighboring subunits.See also Video S2.

Interaction of Pore-Lining Residues with Ions in CLDN10a and CLDN10b Channels
In order to compare the pore-lining surfaces of CLDN10a and CLDN10b channels, the electrostatic surface potential was calculated by solving the Poisson-Boltzmann (PB) equation using the PBEQ Solver (https://charmm-gui.org/?doc=input/pbeqsolver,URL accessed on 30 January 2024 [41]) (Figure 7A-D).Fitting to the above-mentioned opposite net charge of the two pores, the electrostatic potentials differed strongly: Mostly positive for CLDN10a but mostly negative for CLDN10b, especially in the inner half of the pore (Figure 7B,D).The most relevant residues lining the pore (~60 Å in length) are shown in (Figure 6C,D).The pore center is formed by four N54 and four H60 residues belonging to subunits of the tetrameric core barrel of the pore.In contrast, in the CLDN10b model, the pore center was formed by the D56 (negative) and N62 at the corresponding positions [15].Next to these residues, the positively charged R62 and K64 are located in CLDN10a and-10b, respectively.However, close to R62, another positively charged R33 is located in CLDN10a, whereas the negatively charged D36 is located close to K64 in CLDN10b [15].Of note, R33/D36 residues are contributed by the interlocked β1β2 loop of neighboring subunits that are not part of the tetrameric core barrel (octameric-interlocked-barrels architecture [15]).In addition, further, in the direction of the pore periphery in the ECH region, K69 is present in CLDN10a instead of the non-polar A71 in CLDN10b.Both residues, K69 and R33, were capable of neutralizing at least transiently the negatively charged E151/E153 common for CLDN10a and -10b (Figure 3G-J).Furthermore, CLDN10a with Q45, N50, S58, H64, contains more additional polar residues than CLDN10b with the corresponding A47, A52, V60, F66 residues.With these mentioned and other residues (mostly identical for CLDN10a and -10b: K29/31, A71/D73 E143/145, D146/148, K139/141 (K153/K155, E155/157 mainly shielded), the pore has a net charge of +8 for CLDN10a and -12 for CLDN10b.Consequently, the CLDN10a pore was largely filled with Cl − ions, whereas that of CLDN10b with Na + ions (Figures 1 and 6).

Interaction of Pore-Lining Residues with Ions in CLDN10a and CLDN10b Channels
In order to compare the pore-lining surfaces of CLDN10a and CLDN10b channels, the electrostatic surface potential was calculated by solving the Poisson-Boltzmann (PB) equation using the PBEQ Solver (https://charmm-gui.org/?doc=input/pbeqsolver,URL accessed on 30 January 2024 [41]) (Figure 7A-D).Fitting to the above-mentioned opposite net charge of the two pores, the electrostatic potentials differed strongly: Mostly positive for CLDN10a but mostly negative for CLDN10b, especially in the inner half of the pore (Figure 7B,D).
The pore-lining residues were further analyzed by measuring the contact time of relevant residues (see Section 2.3) of CLDN10a and CLDN10b with Na + and Cl − ions.In CLDN10b, the central D56 was predominantly in contact with Na + ions (>80% of the time), D36 (>40% of the time), and E153 (>25% of the time) also showed frequent contact (Figure 7E).In contrast, contacts with Cl − were very rare.In comparison, for CLDN10a, the opposite was obtained: Frequent contacts with Cl − but very rare with Na + .However, for CLDN10a, no residue had a normalized contact time with Cl − ions >45%, and the contacts were spread over more residues than for CLDN10b: The central H60 and N54 with ~34% and ~25% normalized contact time, respectively, and R62, R33, and K69 with ~42%, ~44% and ~23% normalized contact time, respectively.For both channels, charged residues in the inner part of the pore were frequently in contact with oppositely charged ions, whereas charge residues close to the entrance (K139 for CLDN10b, D73 for CLDN10a) had very infrequent contact with the ions (Figure 7E).
Finally, we investigated the charge selectivity of CLDN10 channels by applying an external voltage gradient (Figure 8).The low-resistance pathway on both sides of the dodecamer hindered measurement of the ionic current and calculation of the conductance as determined experimentally [35,39].However, as a quantitative estimate of the ion movement through the channel and as an approximation of the charge selectivity of the pore, total ion displacement between the two ends of the channel over the simulation time was calculated.For CLDN10b, the total displacement of the Na + ions depended linearly on the voltage and was between 578 Å in one direction for −1.4 V and 457 Å in the opposite direction for 1.4 V.In contrast, the total displacement of Cl − ions was 0 Å for all voltages, clearly showing the cation selectivity of the channel (Figure 8A).We also analyzed the CLDN10b mutant K64M that was suggested to enhance Na + ion interaction [15].For CLDN10b-K64M, the slope of the linear fit of the relationship between the total ion displacement and voltage was slightly higher than for CLDN10b-wt, suggesting a slightly higher Na + conductance for the mutant, as expected.For CLDN10a, the total ion displacement of Cl − ions was 544 Å for −1.4 V and 536 Å for 1.4 V in opposite directions.The total ion displacement of Na + ions was −255 Å for −1.4 V and 0 Å for 1.4 V.The absolute value of the slope of the linear fit for the Cl − total displacement was higher than that for the Na + total displacement, indicating the anion preference of the CLDN10a channel (Figure 8B).
for CLDN10a, no residue had a normalized contact time with Cl ions >45%, and the contacts were spread over more residues than for CLDN10b: The central H60 and N54 with ~34% and ~25% normalized contact time, respectively, and R62, R33, and K69 with ~42%, ~44% and ~23% normalized contact time, respectively.For both channels, charged residues in the inner part of the pore were frequently in contact with oppositely charged ions, whereas charge residues close to the entrance (K139 for CLDN10b, D73 for CLDN10a) had very infrequent contact with the ions (Figure 7E). ) and Na + ions.The contact profile differs strongly between the two channels.The residues (CLDN10a/-10b) are ordered according to their position along the x-axis from the pore center towards the two-pore entrances.For the calculation of normalized contact time, the percentage of time frames in which the respective residue was closer than 4 Å to at least one ion was calculated for each subunit lining the central pore and averaged.
In addition, we determined the mean number of ions in the pore during the simulation time (Figure 8C,D).For CLDN10b much more Na + ions were detected than for CLDN10a.For CLDN10b-K64M, even more Na + ions were present.In contrast, many Cl − ions were detected for CLDN10a but hardly any for CLDN10b and CLDN10b-K64M.Interestingly, the number of Cl − ions for CLDN10a and that of Na + ions in CLDN10b were in a similar range.No clear dependence of the mean ion number and the voltage gradient was observed for CLDN10b and Na + as well as CLDN10a and Cl − .However, with decreasing voltage, the Na + ion number increased for CLDN10a, decreased for CLDN10b-K64M, and the Cl − ion number increased slightly with increasing voltage for CLDN10b and CLDN10b-K64M.These voltage dependencies might reflect model asymmetries that are enhanced by strong voltage gradients.However, in sum, the data obtained with the electric field simulations demonstrated opposite charge-dependence of attraction and conduction of ions for the CLDN10a and -10b channel models, respectively.(A) For CLDN10b, the total displacement of Na + (red dots) but not that of Cl − (red triangles) depended linearly on voltage.For the CLDN10b-K64M mutant, the slope of the regression line was higher than for the wild type, suggesting slightly higher conductance for the mutant.(B) For CLDN10a, the absolute value of the slope of the regression line for Cl − total displacement was higher than for Na + total displacement, indicating the opposite charge selectivity of ion conduction compared to CLDN10b.In (A,B), the negative values of total displacement represent the movement in opposite directions on the same axis.(C) The mean number of Na + ions in the pore for CLDN10b was much higher than for CLDN10a and highest for CLDN10b-K64M.(D) For CLDN10a, much more Cl − ions were in the pore than for CLDN10b and CLDN10b-K64M.The mean number of Cl − ions in the pore for CLDN10a was similar to the Na + ion number for CLDN10b.

Discussion
In this study, we investigated paracellular anion channels formed by CLDN10a using homology modeling and MD simulations.We compared the channels with those formed by CLDN10b simulated previously [15] with respect to interfaces between subunits, pore-lining residues, and ion passage.The results indicate that CLDN10a and -10b channels share a common architecture: Interlocked pore barrels as part of a joined double External electric fields were applied as driving forces, and total ion displacement was calculated over simulation time (50 ns).(A) For CLDN10b, the total displacement of Na + (red dots) but not that of Cl − (red triangles) depended linearly on voltage.For the CLDN10b-K64M mutant, the slope of the regression line was higher than for the wild type, suggesting slightly higher conductance for the mutant.(B) For CLDN10a, the absolute value of the slope of the regression line for Cl − total displacement was higher than for Na + total displacement, indicating the opposite charge selectivity of ion conduction compared to CLDN10b.In (A,B), the negative values of total displacement represent the movement in opposite directions on the same axis.(C) The mean number of Na + ions in the pore for CLDN10b was much higher than for CLDN10a and highest for CLDN10b-K64M.(D) For CLDN10a, much more Cl − ions were in the pore than for CLDN10b and CLDN10b-K64M.The mean number of Cl − ions in the pore for CLDN10a was similar to the Na + ion number for CLDN10b.

Discussion
In this study, we investigated paracellular anion channels formed by CLDN10a using homology modeling and MD simulations.We compared the channels with those formed by CLDN10b simulated previously [15] with respect to interfaces between subunits, pore-lining residues, and ion passage.The results indicate that CLDN10a and -10b channels share a common architecture: Interlocked pore barrels as part of a joined double (=quadruple) rows arrangement (JDR) of claudin subunits within TJ strands.In addition to conserved interaction motifs, some inter-subunit interfaces differ between the two isoforms.Furthermore, pore-lining residues differing between CLDN10a and -10b channels cause opposite charge selectivity.The simulations are very well in line with experimental data and provide, for the first time, information on a complete and prototypic paracellular anion channel on the atomic level.
Tracer flux assays and electrophysiological measurements such as transepithelial resistance (TER) or dilution potentials of cellular monolayers are essential to characterize the functionality of TJs constituted by claudins.In combination with mutagenesis and biochemical studies, they can be used to identify sequence determinants for the formation of TJ strands, barriers, and channels [2,42].However, these methods provide only limited information about the structure and molecular dynamics of claudin polymers, TJ strand formation, and how pore-lining residues aid in the ion permeation process.Interestingly, boosted by the first claudin crystal structure solved by Suzuki et al. [9,10], computational modeling and MD simulations have contributed significantly to the understanding of claudin oligomerization and pore formation.
Although different molecular models of claudin-based TJ strands and ion channels have been reported and refined, their validity and differences between the claudin subtypes are disputed [4,12,15,[18][19][20][21][22][23][25][26][27][28][29]43,44].This is also due to the fact that in contrast to transmembrane channels and other cell-junctional structures, no experimental structure has been resolved for TJ-like claudin oligomers so far [4] (https://www.rcsb.org/URL accessed on 10 January 2024).Comparison of the channel-forming CLDN15, CLDN10b, and CLDN10a in this and a previous study [15] by MD simulations support the idea that at least these three claudins form channels and strands according to the JDR arrangement suggested by Suzuki et al. [10].The linear-cis and face-to-face-cis interfaces of this arrangement model are supported by several experimental structure-functions studies on CLDN15, -10b, -5, -3, -2 and -1 [8,[10][11][12]19,24,45,46].Regarding the trans-interfaces, critical contributions by the β1β2 loop of ECS1 and the ECS2 have been proposed, but the interaction patterns were largely unclear [8,10].Modeling and MD simulation studies suggested two different conformational variants by which the β1β2 loop participates in cis-/trans-oligomerization of classic claudins and in conjunction of adjacent tetrameric subunit repeats within strands.In the first variant, the β1β2 loop shows a rather straight extension and orientation towards a subunit in the opposing membrane.This results in a tetrameric-locked-barrel (4LB) that forms a pore in the case of channel-forming claudins and is nearly completely lined by residues from only four subunits [12,15].To a certain extent, similar tetrameric pore variants have been reported earlier [25,28].In the second variant, the β1β2 loop is oriented flatter towards a subunit in the same membrane (Figure 1D).This results in interlocked-barrels (IB) and pores (in the case of channel-forming claudins) that are lined, especially in the inner region, by eight subunits [15,26,27] (Figure 1C,G).Comparison of both variants for CLDN10b and CLDN15 by MD simulations supported the IB model [15].
Consequently, in this study, we simulated IB models for the comparison of CLDN10a and CLDN10b channels.Key results for CLDN10a model variant IB-2, which was chosen as the reference model, and the previously generated CLDN10b IB model [15] are summarized in Table 1.The equilibrated CLDN10a dodecamers showed an overall architecture similar to that of CLDN10b.This was expected since a CLDN10b template was used to generate the CLDN10a homology dodecamer models.More importantly, the CLDN10a model variants showed similar stability as the CLDN10b models: RMSD of dodecamer protein backbone was similar (Figure S4), contact map pattern (Figures 2 and 4), as well as interaction parameter for face-to-face, linear-cis, ECS2-ECS2-trans (Figure 3) and hydrophobic cluster of β1β2 loop tips (Figure 5), were in most cases slightly higher for CLDN10a models than for CLDN10b models [15], however still in similar ranges.All interfaces and an open pore conformation were preserved during the simulations (Table 1).Thus, the results support the IB model for CLDN10a and suggest that at least CLDN15, -10b, and -10a share a similar overall architecture: Adjacent, sidewise-unsealed tetrameric pore scaffolds that are interlocked via β1β2 loops.These loops mediate, together with ECS2, cisand trans-interaction between subunits of adjacent tetramers.The clustering of the hydrophobic β1β2 loop tips of four subunits was proposed to be a common driving force for the oligomerization of classic claudins for which the hydrophobicity is conserved [4,7,13].
Strikingly, the ion conduction pathway differed strongly between CLDN10a and -10b in the inner half of the pore (Figures 1 and 6).While the pore center in CLDN10a models is formed by N54 and H60 residues, it is formed by D56 and N62 residues at the corresponding positions in CLDN10b models.Close to the center, both claudins contain a different positively charged residue (R62/K64).Next to this position, CLDN10a contains another positively charged R33, whereas in CLDN10b, the negatively charged D36 is located.Importantly, these residues are contributed by the interlocked β1β2 loop of adjacent subunits that are not part of the tetrameric core barrel.Further to the pore periphery, the charged K69 is present in CLDN10a instead of the non-polar A71 in CLDN10b.CLDN10a also contains more additional polar pore-lining residues than CLDN10b.In sum, this leads to an opposite electrostatic potential on the pore-lining surface for the two different claudin channels (Figure 7A-D).Consequently, the pore-lining residues of CLDN10a attracted and interacted strongly with Cl − ions, whereas that of CLDN10b did so with Na + ions (Figures 7E and 8C,D).This resulted in anion transport through CLDN10a channels in contrast to cation transport through CLDN10b channels (Figure 8).Thus, in total, the simulations fit very well to the charge-selectivity of the two claudin channels and contributions of ECS1 residues, in particular of R33 and R62 in CLDN10a, that were demonstrated by electrophysiological measurements [35,36,38,39].
As mentioned above, the equilibrated CLDN10a and CLDN10b models were further used to study the ion permeation.However, it is important to mention that reproducing the physiological ion transport through claudin channels using computational models is complicated due to various possible sources of error [51].During an all-atom MD simulation of a claudin pore model, movement of only a few ions through the pore could be observed because of the limited simulation time.Moreover, since these paracellular pores do not run through a membrane barrier, it is very challenging to calculate a net flux in such simulations.Different previous computational studies on claudin channel models approached this problem in different ways.Alberini et al. [43] and Irudayanathan et al. [18,52] observed the traverse of different ions through claudin channel models by calculating the potential of mean force (PMF) profiles using Umbrella sampling simulations.Samanta et al. [26] calculated the total ionic currents of cations through CLDN15 channels after applying a constant electric field corresponding to a transepithelial potential during the simulation.We followed a similar methodology in this study; however, there are some important limitations.One important issue is the use of conventional periodic boundary conditions (PBC), which results in a continuous bulk solution; hence, the ions are free to diffuse through boundaries.Samanta et al. countered this problem by allowing the claudin strands to continue across the PBCs, as a non-continuous strand (as in our system setups) will result in a low-resistance pathway on the two ends of the periodic box, which results in leakage of ions.To address this, we concentrated only on the ion movement through the middle pore and not the whole system by considering an imaginary cuboid through the pore.The ion movement events occurring in the cuboid were observed and quantified using displacement calculations.Although this approach hinders the calculation of current and conductance through the channels across an epithelium as measured experimentally [35,39], it allows a relative comparison of ion movements through the different claudin channels.
In total, this MD simulation study provides the first atomic model of a complete and prototypic paracellular anion channel and suggests that CLDN10a anion and -10b cation channels share the same overall architecture but differ mainly in their pore-lining residues.This is of high interest since nearly all transcellular anion channels differ strongly from cation channels in their architecture [53].However, one has to keep in mind that claudinbased paracellular channels differ strongly in structure from transmembrane channels since the ion pore runs parallel to instead of through the membrane.A similar overall architecture for paracellular anion and cation channels also fits the fact that both are-in contrast to transmembrane channels-embedded in protein polymers that form a functional barrier against larger solutes.
Certainly, the first structural model reported here cannot reflect all structural and mechanistic details of native CLDN10a channels.On the one hand, more structure-function studies employing mutagenesis and electrophysiology have to be performed to identify more sequence determinants for defined channel properties.On the other hand, more detailed MD simulations that take the in vitro data into account have to be performed.Varying details of regional conformations in the starting models, the lipid composition of the membrane, ions, force fields, release of remaining constraint (C61), number of subunits, and longer simulation times can be used for structural refinement in further studies.Thereby, for instance, it is expected to reduce channel asymmetries and, in turn, increase the accuracy of the simulations.In addition, free energy calculation can be used to analyze selectivity determinants in more detail [18,21,23,43].Furthermore, progress in claudin complex structures resolved by cryo-electron microscopy [53] is expected to provide additional substantial insights into the structure of claudin channels and strands.
In sum, we generated a structural model for the prototypic anion-selective CLDN10a channel.Molecular dynamics simulations of CLDN10a or CLDN10b dodecamers indicate that both CLDN10 isoforms share the same channel and strand architecture: Sidewise unsealed tetrameric pore scaffolds, interlocked with adjacent pores via β1β2 loop of ECS1, leading to strands with claudin subunits arranged in double rows in one and in four joined rows in two membranes.Many interaction modes are suggested to be conserved among CLDN10a, -10b, and also -15.However, pore-lining residues differing between CLDN10a/-10b (i.e., R33/I35, A34/D36, K69/A71, N54/D56, H60/N62, R62/K64) result in opposite charge selectivity of the two channels.The findings are consistent with electrophysiological studies and thus provide novel mechanistic information about epithelial transport.

Modeling of CLDN10a IB Dodecamer Model
A homology model of human CLDN10a (Uniprot P78369-2, M1-N184) was developed using the 'Build Homology model' module in Maestro BioLuminate.A well-equilibrated CLDN10b dodecamer subunit [15] was used as the modeling template.In addition, the conformation of the ECH region (66-71) was grafted from a CLDN10a model predicted using ColabFold (https://colab.research.google.com/github/sokrypton/ColabFold/blob/main/AlphaFold2.ipynb,URL accessed on 5 April 2023) [54,55].An initial CLDN10a dodecamer model was modeled similarly to the previously published CLDN10b dodecamer models by replicating the CLDN10a homology model and aligning it with the CLDN10b dodecamer [15].The CLDN10a dodecamer model consisted of two trans-interacting hexamers, resulting in a triple-pore model.A series of minimizations based on different experimental-and hypothesis-derived distance constraints [4,12,15,38,39] were performed using the 'Macromodel minimization' tool available in Maestro BioLuminate.The Polak-Ribiere Conjugate Gradient (PRCG) method that uses the Polak-Ribiere first derivative method was used [56].The RMS gradient of the energy with respect to the coordinates in kJ mol −1 Å −1 was selected as the convergence criterion.After minimization, a series of molecular dynamics (MD) simulations in a water-solvated environment was followed to pre-equilibrate the model, considering the same distance constraints.The following considerations based on conservation of (a) sequence motifs and (b) strand formation among classic claudins [4,12,15,38,39] were followed during model equilibration: i.
The CLDN10a-specific, positively charged residues R33 and K69 are oriented towards the anion-selective pore.

MD Simulations of CLDN10a IB Dodecamer Models
'Desmond Molecular Dynamics' module (Desmond Molecular Dynamics System, D. E. Shaw Research, New York, NY, USA, 2022.Maestro-Desmond Interoperability Tools, Schrödinger, New York, NY, USA, 2022; [15,57]) present in Maestro BioLuminate was used to perform the MD simulations of the CLDN10a dodecamer.The pre-equilibrated CLDN10a IB dodecamer model was embedded in two 1-palmitoyl-2-oleoyl-sn-glycero-3phosphocholine (POPC) lipid bilayers.An equilibrated membrane from the CLDN10b 8IB (short: IB) model [15] was used to follow a similar positioning of the dodecamer in the membranes.The membrane was further refined to remove lipids clashing with protein and to add missing lipids to avoid water inclusion.The double membrane-dodecamer complex was put into a simulation box of size 165, 165, and 150 Å in the x-, y-and z-axis, respectively.The orthorhombic box was solvated with TIP3P water molecules [58], physiological salt of 0.15 M Na + Cl −, and charge-neutralizing Cl − ions.OPLS4 force field [59] was applied during the simulations performed in the NPT ensemble at 310 K and 1.01325 bar.Nosé-Hoover chain method [60] and Martyna-Tobias-Klein method [61] were used as thermostats and as barostats, respectively.The short-range cutoff method with a cutoff radius of 9.0 Å was used to calculate Coulombic interactions.Bonded interactions were integrated using the RESPA algorithm at 2.0 fs timestep.The system was first relaxed into a local energy minimum through a 'Desmond minimization' available in Maestro BioLuminate, using a Brownian motion simulation run for 100 ps, followed by the NPT BioLuminate default relaxation protocol with added constraints on protein except for H atoms. Subsequently, a series of equilibration steps were performed for 160 ns by gradually releasing the constraints.In the beginning, a force constant of 10 kcal mol −1 Å −2 was applied to the whole protein.Afterward, first side chains, then backbone atoms, stepwise for loops, β-sheets, and helices were gradually released over time by decreasing the respective force constants stepwise from 5 to 0 kcal mol −1 Å −2 .
Additionally, one of the predecessors of the above-mentioned CLDN10a IB-1 model was refined further to increase the interface-wise symmetry of the model by adjusting some interfacial residues (such as S35, S36, L56, F68, P147, Q152) using 'Maestro builder'.The refined model IB-2 was embedded in a double membrane, minimized, stepwise equilibrated, and simulated similarly to the IB-1 model.For analysis of the IB-2 model, four variant production runs were performed that differed in the following constraints: i.
Modeling and simulation steps that lead to the above-mentioned model variants are summarized in Figure S3.

Analysis of the MD Trajectories
We extracted relevant data, including RMSD, interactions like hydrogen bonds and salt bridges, and distance between residues from simulation trajectories using the internal tools available in BioLuminate.The extracted data were compiled and analyzed using Microsoft excel (Microsoft Corporation, 2018, Microsoft Excel, Redmond, WA, USA).The analysis of Solvent-Accessible Surface Area (SASA) was performed using Visual Molecular Dynamics [62] (VMD, version 1.9.4a55, release 2021-10, https://www.ks.uiuc.edu/Research/vmd,IL, USA).
MDAnalysis (version 2.3.0)[63,64] scripts written using Python3 [65] were developed to perform the analysis and visualization of pore dimension using the HOLE program [66,67] and residue-residue interaction analysis using contact maps.More information on the scripts can be found in our previous publication [15].

Applied Electric Field Simulations
We performed simulations by applying a potential bias in the axis parallel to the membrane surface and through the paracellular pores.One of the well-equilibrated structures from the production runs of CLDN10a (IB-2) and CLDN10b (8IBli) [15] simulations was selected.The protein dodecamer structures, along with the membranes, ions, and water, were converted into a new simulation system, energy minimized, and then equilibrated for 5 ns.During equilibration, the protein and the head group of lipid molecules were constrained with a force constant of 1 kcal mol −1 Å −2 .Using the 'e_bias' plugin available in the Desmond simulation platform, simulations with applied electric fields along the membrane surface with different potentials like 1.4 V, 0.8 V, 0.4 V, −0.4 V, −0.8 V, and −1.4 V were simulated.During the simulations, the extracellular segments were kept free, and only the backbone atoms of the transmembrane helices were constrained with a force constant of 1 kcal mol −1 Å −2 to avoid their movement in the direction of the current due to the sudden application of the potential.The production runs were simulated for 50 ns in the NVT ensemble.

Calculation of Ion Total Displacement
As a quantitative measure for the ion permeation, using MDAnalysis, the total displacement of Na + and Cl − ions through the pore was calculated using the coordinate information from the trajectories.The displacement calculation was performed for only the ions that cross a bottleneck in the center of the middle pores of CLDN10a and CLDN10b dodecamers.For the bottlenecks, the tetrads of H60 in CLDN10a and D56 in CLDN10b were considered, respectively, as these residues from four different chains come together at the center of the pore pathway.So, any ion that passes the center of mass of these four residues was considered to be passing the pore.However, a complication in choosing the ions traveling and passing only through the pore and not through any other side arose.For this purpose, we set boundaries in x, y, and z-axes to consider an imaginary cuboid through the pore, and only the ions that pass the pore within the cuboid were selected for the displacement calculation.The Cartesian coordinates of the C α atom of K139 (K141 in CLDN10b) from the four chains that are located at the entrance of the middle pore on either side were used as the boundaries of the cuboid (Figure S6).The displacements of the filtered ions were calculated by extracting the position information of these ions along the

Figure 1 .
Figure 1.Comparison of CLDN10b and CLDN10a dodecamer (three pore) models.(A-C) Snapshots of the production run at 100 ns for the CLDN10b model (CLDN10b_8IBno, taken and modified from [15]).(A) Front view; (B) turned view on the middle pore, claudin subunits (chains) shown as a colored cartoon, POPC lipids as lines, and phosphate headgroups as spheres.(E-G) Snapshot of the production run at 100 ns for CLDN10a model (IB-2): (E) front view;(F) turned view on the middle pore.Protein chains were embedded well in membranes, and the overall arrangement, including trapped lipids (orange) between the two claudin rows, was similar to that of CLDN10b.In the CLDN10a pore, mainly Cl − ions (green spheres) were present, whereas in that of CLDN10b, mainly Na + ions (magenta spheres) were present, fitting to opposite charge selectivity of the two different channels.(D) Schema, showing the contact of β1β2 loops (red and green) between trans-interacting subunits (other claudin parts gray), highlighting the kinked/flat orientation of the loops towards subunits in the same membrane.This results in interlocked loops of neighboring pores (interlocked barrels (IB) arrangement).CDLN10a and -10b models showed similar IB arrangements.(C,G) Close-up of dodecamer centers.On both sides of the middle pore, hydrophobic clusters were formed similarly for CLDN10a and -10b by V37/I38 (CLDN10a) and V39/I40 (CLDN10b) residues (shown as surfaces) of trans-and cis-interacting β1β2 loop tips from four subunits (blue, green, red, and black).In contrast, hydrophilic pore-lining residues (sticks) differed strongly between CLDN10a and -10b.See also Video S1.

Figure 1 .
Figure 1.Comparison of CLDN10b and CLDN10a dodecamer (three pore) models.(A-C) Snapshots of the production run at 100 ns for the CLDN10b model (CLDN10b_8IBno, taken and modified from [15]).(A) Front view; (B) turned view on the middle pore, claudin subunits (chains) shown as a colored cartoon, POPC lipids as lines, and phosphate headgroups as spheres.(E-G) Snapshot of the production run at 100 ns for CLDN10a model (IB-2): (E) front view; (F) turned view on the middle pore.Protein chains were embedded well in membranes, and the overall arrangement, including trapped lipids (orange) between the two claudin rows, was similar to that of CLDN10b.In the CLDN10a pore, mainly Cl − ions (green spheres) were present, whereas in that of CLDN10b, mainly Na + ions (magenta spheres) were present, fitting to opposite charge selectivity of the two different channels.(D) Schema, showing the contact of β1β2 loops (red and green) between trans-interacting subunits (other claudin parts gray), highlighting the kinked/flat orientation of the loops towards subunits in the same membrane.This results in interlocked loops of neighboring pores (interlocked barrels (IB) arrangement).CDLN10a and -10b models showed similar IB arrangements.(C,G) Closeup of dodecamer centers.On both sides of the middle pore, hydrophobic clusters were formed similarly for CLDN10a and -10b by V37/I38 (CLDN10a) and V39/I40 (CLDN10b) residues (shown as surfaces) of transand cis-interacting β1β2 loop tips from four subunits (blue, green, red, and black).In contrast, hydrophilic pore-lining residues (sticks) differed strongly between CLDN10a and -10b.See also Video S1.

Figure 2 .
Figure 2. Inter-subunit residue contact maps for CLDN10b and CLDN10a models.(A) Overview of different inter-subunit interfaces.Middle: Front view on CLDN10a dodecamer with subunits as colored cartoons.Left: Turned view on ECS1+2 of two bottom subunits (dashed blue box in the middle) highlighting face-to-face cis interface with H-bonds between β4 strands.Right: Close-up of the dashed black box in the middle highlighting the linear-cis interface between the ECH region (blue) and ECS2 pocket (green) and trans-interface between two ECS2 turn regions (green, red).Relevant residues are shown as sticks.Far right: Comparison of CLDN10a and CLDN10b ECH regions.Overlay of several subunits showing different orientations of F68/I67 in CLDN10a and of corresponding L70/M69 in CLDN10b model.(B-F) Mean distances (closest atoms) between the numbered residues of protein region pairs reflecting the different interface types.For each interface type, the multiple individual interfaces in the dodecamer were averaged.Shown are interaction-relevant parts of contact maps categorized according to interface types: (B) face-to-face-cis, (C) linear-cis, (D) ECS2-ECS2-trans, and (E,F) β1β2 loop-ECS2.In headlines, corresponding residue numbers of CLDN10a are given in brackets.Residue contact pairs are boxed in white if similar, in red if shifted, and in black when different between CLDN10a and -10b.The contact maps provide interface fingerprints for comparison of the different claudin models.CLDN10b data taken from previous simulations [15].

Figure 2 .
Figure 2. Inter-subunit residue contact maps for CLDN10b and CLDN10a models.(A) Overview of different inter-subunit interfaces.Middle: Front view on CLDN10a dodecamer with subunits as colored cartoons.Left: Turned view on ECS1+2 of two bottom subunits (dashed blue box in the middle) highlighting face-to-face cis interface with H-bonds between β4 strands.Right: Close-up of the dashed black box in the middle highlighting the linear-cis interface between the ECH region (blue) and ECS2 pocket (green) and trans-interface between two ECS2 turn regions (green, red).Relevant residues are shown as sticks.Far right: Comparison of CLDN10a and CLDN10b ECH regions.Overlay of several subunits showing different orientations of F68/I67 in CLDN10a and of corresponding L70/M69 in CLDN10b model.(B-F) Mean distances (closest atoms) between the numbered residues of protein region pairs reflecting the different interface types.For each interface type, the multiple individual interfaces in the dodecamer were averaged.Shown are interaction-relevant parts of contact maps categorized according to interface types: (B) face-to-face-cis, (C) linear-cis, (D) ECS2-ECS2-trans, and (E,F) β1β2 loop-ECS2.In headlines, corresponding residue numbers of CLDN10a are given in brackets.Residue contact pairs are boxed in white if similar, in red if shifted, and in black when different between CLDN10a and -10b.The contact maps provide interface fingerprints for comparison of the different claudin models.CLDN10b data taken from previous simulations [15].

Figure 4 .
Figure 4. Contact maps for interfaces between different pairs of β1β2 and β3β4 loops in CLDN10b and CLDN10a models.(A) Overview.Left: Front view on CLDN10a dodecamer.Middle: A close-up of the dashed black box on the left shows a central cluster of different loops of neighboring subunits (different colors) in contact.Right: Turned view on this β1β2/β3β4 loop cluster showing another perspective on loop contacts.(B-D) Mean distances (closest atoms) between the numbered residues of β1β2 loops (CLDN10a: 30-45) and β3β4 loops (CLDN10a: 52-58) of three different subunit pairs.For the definition of pairs, see headlines and subunit numbering according to (A).The multiple individual interfaces in dodecamer were averaged.Residue contact pairs are boxed in white if similar, in red if shifted, and in black when different between CLDN10a and -10b.The contact maps provide interface fingerprints for comparison of the different claudin models.CLDN10b data taken from previous simulations [15].

Figure 4 .
Figure 4. Contact maps for interfaces between different pairs of β1β2 and β3β4 loops in CLDN10b and CLDN10a models.(A) Overview.Left: Front view on CLDN10a dodecamer.Middle: A close-up of the dashed black box on the left shows a central cluster of different loops of neighboring subunits (different colors) in contact.Right: Turned view on this β1β2/β3β4 loop cluster showing another perspective on loop contacts.(B-D) Mean distances (closest atoms) between the numbered residues of β1β2 loops (CLDN10a: 30-45) and β3β4 loops (CLDN10a: 52-58) of three different subunit pairs.For the definition of pairs, see headlines and subunit numbering according to (A).The multiple individual interfaces in dodecamer were averaged.Residue contact pairs are boxed in white if similar, in red if shifted, and in black when different between CLDN10a and -10b.The contact maps provide interface fingerprints for comparison of the different claudin models.CLDN10b data taken from previous simulations[15].

Figure 6 .
Figure 6.Ion permeation pathway of pore in CLDN10a models.(A,C,D) Snapshots at 100 ns for IB-2 model.Eight subunits lining the middle pore are shown as colored cartoons; the permeation pathway determined by the HOLE program is shown as a transparent gray surface, Na + as red, and Cl − as green spheres.(A) Overview: The pore axis runs from left to right.(B) Pore diameter profiles for CLDN10a IB-2, IB-2+lic, IB-2+hc and IB-2+ohc models.The pore profile of IB-2 is represented by a bold black line, with the corresponding SD values given as gray bars in addition.Mean diameter of last 50 ns of simulation along the pore axis.Pore pathway and diameter detection by HOLE.The position of the most relevant residues along the pore axis is indicated.(C,D) Close-up of permeation pathway in two different orientations.Most-relevant pore-lining residues are shown as sticks.R33 residues belong not to the tetrameric pore scaffold but to neighboring subunits.See also Video S2.

Figure 6 .
Figure 6.Ion permeation pathway of pore in CLDN10a models.(A,C,D) Snapshots at 100 ns for IB-2 model.Eight subunits lining the middle pore are shown as colored cartoons; the permeation pathway determined by the HOLE program is shown as a transparent gray surface, Na + as red, and Cl − as green spheres.(A) Overview: The pore axis runs from left to right.(B) Pore diameter profiles for CLDN10a IB-2, IB-2+lic, IB-2+hc and IB-2+ohc models.The pore profile of IB-2 is represented by a bold black line, with the corresponding SD values given as gray bars in addition.Mean diameter of last 50 ns of simulation along the pore axis.Pore pathway and diameter detection by HOLE.The position of the most relevant residues along the pore axis is indicated.(C,D) Close-up of permeation pathway in two different orientations.Most-relevant pore-lining residues are shown as sticks.R33 residues belong not to the tetrameric pore scaffold but to neighboring subunits.See also Video S2.

Figure 7 .
Figure 7. Electrostatic surface potential (ESP) map of CLDN10a (A,B) and CLDN10b (C,D) channel models.(A,C) Electrostatic potential mapped on dodecamer protein surface of 100 ns snapshots of CLDN10a IB-2 and CLDN10b [15] simulations.Front view along y and z axes.(B,D) Clipped cross-section along x and z axes at the middle position of dodecamers (dashed line in (A,C)).The internal protein region is shown in black.Electrostatic surface potential calculated with Poisson-Boltzmann equation Solver (PBEQ-Solver) is shown color-coded (−10 to +10 kCal/mol.e).ESP of the pore-lining surface differs strongly between CLDN10a and -10b.(E) Contact time of relevant pore-lining residues with Cl −) and Na + ions.The contact profile differs strongly between the two channels.The residues (CLDN10a/-10b) are ordered according to their position along the x-axis from the pore center towards the two-pore entrances.For the calculation of normalized contact time, the percentage of time frames in which the respective residue was closer than 4 Å to at least one ion was calculated for each subunit lining the central pore and averaged.

Figure 8 .
Figure8.CLDN10a and -10b channel models show opposite charge selectivity.External electric fields were applied as driving forces, and total ion displacement was calculated over simulation time (50 ns).(A) For CLDN10b, the total displacement of Na + (red dots) but not that of Cl − (red triangles) depended linearly on voltage.For the CLDN10b-K64M mutant, the slope of the regression line was higher than for the wild type, suggesting slightly higher conductance for the mutant.(B) For CLDN10a, the absolute value of the slope of the regression line for Cl − total displacement was higher than for Na + total displacement, indicating the opposite charge selectivity of ion conduction compared to CLDN10b.In (A,B), the negative values of total displacement represent the movement in opposite directions on the same axis.(C) The mean number of Na + ions in the pore for CLDN10b was much higher than for CLDN10a and highest for CLDN10b-K64M.(D) For CLDN10a, much more Cl − ions were in the pore than for CLDN10b and CLDN10b-K64M.The mean number of Cl − ions in the pore for CLDN10a was similar to the Na + ion number for CLDN10b.

Figure 8 .
Figure8.CLDN10a and -10b channel models show opposite charge selectivity.External electric fields were applied as driving forces, and total ion displacement was calculated over simulation time (50 ns).(A) For CLDN10b, the total displacement of Na + (red dots) but not that of Cl − (red triangles) depended linearly on voltage.For the CLDN10b-K64M mutant, the slope of the regression line was higher than for the wild type, suggesting slightly higher conductance for the mutant.(B) For CLDN10a, the absolute value of the slope of the regression line for Cl − total displacement was higher than for Na + total displacement, indicating the opposite charge selectivity of ion conduction compared to CLDN10b.In (A,B), the negative values of total displacement represent the movement in opposite directions on the same axis.(C) The mean number of Na + ions in the pore for CLDN10b was much higher than for CLDN10a and highest for CLDN10b-K64M.(D) For CLDN10a, much more Cl − ions were in the pore than for CLDN10b and CLDN10b-K64M.The mean number of Cl − ions in the pore for CLDN10a was similar to the Na + ion number for CLDN10b.

Table 1 .
Comparison of key parameters of CLDN10a and CLDN10b channel models.