Structural dynamics of tight junctions modulate the properties of the epithelial barrier

Tight junctions are dynamic structures that are crucial in establishing the diffusion and electrical barrier of epithelial monolayers. Dysfunctions in the tight junctions can impede this barrier function and lead to many pathological conditions. Unfortunately, detailed understanding of the non-specific permeation pathway through the tight junctions, the so-called leak pathway, is lacking. We created computational models of the leak pathway to describe the two main barrier measures, molecular permeability and transepithelial electric resistance while using common structural dynamics. Our results showed that the proposed alternatives for the leak pathway, the bicellular strand opening dynamics and the tricellular pores, contribute together with distinct degrees, depending on the epithelium. The models can also capture changes in the tight junction barrier caused by changes in tight junction protein composition. In addition, we observed that the molecular permeability was markedly more sensitive to changes in the tight junction structure and strand dynamics compared with transepithelial electric resistance. The results highlight that our model creates a good methodological framework to integrate knowledge on the tight junction structure as well as to provide insights and tools to advance tight junction research.


Introduction
Epithelial cell monolayers cover body surfaces and line different organs. These tissues separate the underlying organs from their surroundings by creating tight barriers, and cell-cell junctions play a crucial role in this process. The most significant components for the barrier function are the tight junctions (TJs). These dynamic structures bring the membranes of adjacent cells into close contact, and thus seal the paracellular space between them. Due to their important role in the epithelial function, it is not surprising that there are several diseases, such as inflammatory bowel disease and celiac disease, which are linked to dysfunctions in TJ proteins or in the TJ complexes themselves [1,2]. In these pathological conditions, the epithelium usually becomes leaky [3], and thus rendering it unfit for its task. In the present work, we a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 the TER increases. The permeability of noncharged molecules has been approached by studying how small tracer molecules, such as polyethylene glycols (PEGs) [31][32][33][34][35][36], dextrans [16,37,38], and mannitol [33,39], traverse the epithelial barrier. Especially different sizes of PEGs have been an invaluable tool when studying how molecular size (<1 kDa) affects the TJ permeation [31][32][33][34][35][36], since they have been shown to permeate paracellularly [40,41].
Molecules and ions are hypothesized to permeate through the TJs via two routes: pore and leak pathways [31,33,42,43]. Although these pathways were defined based on the permeability of noncharged molecules [31,33], they can be extended to the context of the TER. The pore pathway is a high-throughput pathway for molecules with a radius smaller than 0.4 nm through pores formed by some claudins [31][32][33]. The ion permeability of these pores depends on the ion-specificity of the claudins forming them [44]. The leak pathway is a nonspecific, low-throughput pathway with an assumed size limit of over 6 nm [20,31,33,45]. The origin of this pathway is debated, the two main candidates being 1) the pores in the tTJ central tubes and 2) transient bTJ strand breaks [16,20,22,33,[45][46][47][48]. The tTJ pores have been shown to enable macromolecule diffusion: Krug et al. [16] showed that 3-kDa dextran passes through the epithelium via the tTJ pores and suggested that they could form the leak pathway. They further showed that changes in the expression of tricellulin, a TJ protein found especially in the tTJs, affected the permeability of macromolecules, but only minor changes in the permeability of small (<1 kDa) molecules or in TER were seen [16].
The bTJ strands have been shown to remodel constantly at the protein-level [43,49] as well as at the strand-level in transfected fibroblasts by constant strand breaking and sealing events [46,50]. The bTJ strand breaks have also been observed in freeze fracture electron microscope images [16,[51][52][53]. Furthermore, Van Itallie and coworkers recently showed that the connection between ZO-1 and claudins stabilizes the movement of the strands while, interestingly, not affecting the break frequency [50]. In addition, ZO-1 knockdown and double ZO-1/2 knockdown have been demonstrated to increase the leak pathway permeability while only slightly affecting the TER [34,54]. Based on these findings, it has been suggested that structural dynamics in the strands could enable a step-by-step passage for molecules too large to pass through claudin pores. However, the slow dynamics in the strand structure would not be visible in the almost-instantaneous TER measurement [22,33,34,46,48]. Moreover, Liang & Weber [20] suggested that these two pathways are not mutually exclusive.
In recent years, experimental research on the TJ structure and function has been abundant and many important advances have been made [14,18,[55][56][57][58]. However, the small size of the TJs and the molecules passing through them, together with the fast time scales of the events, make TJs and especially their structural dynamics challenging to study experimentally. Computational models complement experimental work and provide an excellent tool to investigate the characteristics of the TJ barrier in more detail. Most of the previous computational models of the TJ barrier [47,[59][60][61][62][63][64] describe a simplified and static TJ structure and do not include any level of structural dynamics. Weber and coworkers have constructed models describing the transient opening behavior of the claudin pores without the strand-level dynamics [65,66]. To the best of our knowledge, the only computational model to include the strand dynamics is a percolation analysis model of the TJ strands as a random resistor network by Washiyama et al. [67]. However, a model that describes the structural TJ dynamics in relation to both the molecular permeability and the TER seems to be lacking.
In this work, we developed computational models for both the molecular permeability and the TER with common structural dynamics to investigate the origin of the leak pathway. We fitted the models to experimental data from two strains of MDCK monolayers and Caco-2 monolayer to parametrize the TJ structure. This, in addition to studying the leak pathway, enabled us to investigate how the different properties of this dynamic system affect the molecular permeability and the TER. With our models, we aim to fill the gap in knowledge between the structural and the functional properties of the TJs.

Modeling framework
The computational models of the dynamic TJ structure for both the molecular permeability and TER use common geometry and dynamic parameters. Both models are divided into the bTJ and the tTJ components. The bTJ geometry is constructed as a two-dimensional structure with a uniformed depth based on the vertical model plane cutting through the TJs in the thin lateral space between the cells as shown in Fig 1B. The real TJ strand network structure is further simplified to an overlapping tile-like ordered network (Figs 1C and 2). Only a section of the strand network is modeled, but the models are averaged for the whole epithelium by long simulation times and by running them multiple times. The main assumptions of the models are as follows: • TJs form the governing barrier against ionic and permeation of noncharged molecules in the epithelium • Leak pathway is formed by both the static tricellular pores and by the bicellular strand dynamics, as suggested by Liang & Weber [20] • tTJ pores are assumed to be similar between epithelia • bTJ strands are impermeable to molecules too large to pass through claudin pores but not to ions • bTJ strands undergo stochastic breaking and resealing events in the scale of seconds to tens of seconds • bTJ strands have homogeneous properties throughout the network • In the time scale of the model, molecular diffusion rate and fluid resistance inside the compartments have no effect on permeability and TER, respectively.
The bTJ molecular permeability model is based on a multi-compartmental approach. Thus, the strand dynamics and the geometry presented as an example structure in Fig 2A are incorporated into rate constants that describe the permeation between the compartments lined by the bTJ strands. The rate constants depend on a stochastic component that describes the strand state as either intact or broken, and whose value changes based on given probabilities. The amount of substance in the basal compartment below the strand network remains constant and the amount of substance in the apical compartment above the strand network is used to calculate the molecular permeability. The model is simulated using a PEG molecule with the mass of 547 Da (calculated radius of 0.51 nm), since it is unable to permeate through the claudin pores and has been used in many experimental permeability studies. The static tTJ pore pathway is combined with the bTJ results afterwards. The bTJ resistance model uses the same geometry and dynamics as the molecular permeability model. However, the system is solved using nodal analysis and Kirchhoff's circuit laws. An example resistor network system is shown in Fig 2B. Current loops form in the resistor network as depicted, and Kirchhoff's loop rule is used for each loop, combined into an equation group, and solved. The dynamics and the geometry shown in Fig 2B are incorporated into the resistances between the compartments. In a similar way to the permeability model, the tTJ resistance is calculated separately and summed with the bTJ simulation results.

Molecular permeability model
The permeabilities of the bTJ and the tTJ components are calculated separately and then combined in the end based on the parallel connection. The bTJ molecular permeability model is constructed using a multi-compartmental approach. The amount of substance in the small bTJ compartments as a function of time is described by an equation of the form where q i (t) is the amount of substance in compartment i at time t and k ji (t) is the time-dependent rate constant for permeation from compartment j to i (s -1 ). Because the model is in 2D, the unit of q i is m -1 instead of unity. This is taken into account when the results are calculated with Eq 8. Since we consider the apical compartment to be considerably larger than the small bTJ compartments, we assume that the concentration in it remains very low. Thus, we can ignore the backflow of molecules from the apical compartment into the bTJ compartments, as shown in Fig 2A. The amount of substance in the basal compartment is assumed to remain constant (dq basal (t)/dt = 0).
Because we assume no molecular permeation through an intact strand, the rate constants for the permeation between the compartments depends only on the strand break permeability coefficient: where l break is the size of the break in the strand, as indicated in Fig 2C (m), r ij (t) is a function describing the state of the strand between compartments i and j, A i is the area of compartment i (m 2 ), and P break is the permeability coefficient of the strand break of indefinite length (m s -1 ). In a 3D model, break area and compartment volume would be used instead of the break size and compartment area, respectively [68]. Function r ij (t) can have a value of either 0 or 1, describing intact and broken states, respectively. The value of r ij (t) can change with given probabilities: If the strand is intact, a break will form with the probability p break (m -1 s -1 ), and if the strand is broken, it will seal with the probability p seal (s -1 ). To obtain the break forming probability for a section of strand between two compartments, p break must be multiplied by the length of that section. This means that longer sections of strand have higher value of p break . Based on the time scale of the dynamics, the possible changes in the strand states were chosen to happen every second.
The initial state of the strands can be either intact or broken. The probability of them being in either state is calculated with Markov chain after infinite time: where p ij,intact and p ij,broken are the probabilities of the strand section between compartments i and j being in intact or broken states initially, respectively, and l ij is the length of the strand section between those compartments (m).
Since the break size is taken into account in Eq 2, P break is calculated relative to the amount of TJs in the entire epithelium, and thus making it dependent on the cell boundary length per area of epithelium: where � TJ is the relative area of the TJs in the epithelium, D 0 is the aqueous diffusion coefficient of the permeating molecule (m 2 s -1 ), H s (r m /w TJ ) is the break slit hindrance factor that depends on the molecular radius r m (m) and on the TJ half-width w TJ (m), and h strand is the bTJ strand height (m) (See Fig 2C for an illustration of the geometrical parameters) [69]. Parameter � TJ is calculated as where l cb is the cell boundary length per area of epithelium (m -1 ). Function H s describes how the break walls affect the diffusion of a molecule of a given size and was derived by Dechadilok & Deen [70] by fitting a polynomial to computational results as where λ = r m /w TJ . With zero initial conditions for the bTJ compartments, there is a so-called lag phase in the beginning of a simulation, especially with systems having more horizontal strands and lower value of p break . During this phase, the increase in the amount of substance in the apical compartment (q apical ) is nonlinear. Since the permeability is calculated from the linear phase in q apical , the simulation is made to enter straight into this phase by setting the initial values of the amount of substance in each of the small bTJ compartments to the equilibrium state during the linear phase. This is done by simulating the permeability model multiple times and for a long time with zero initial conditions to obtain the equilibrium values for each compartment row.
The constant amount of substance for the basal compartment is set to where c basal is the basal compartment concentration (M), N A is Avogadro's constant (6.022 × 10 23 mol -1 ), and A basal is the area of the basal compartment, replacing the volume since the model is in 2D (m 2 ). The system of differential equations described by Eq 1 is solved using Matlab's (Release 2015b, The MathWorks, Inc., Natick, Massachusetts, United States) ode23 ordinary differential equation solver. This solver uses second and third order Runge-Kutta formulas, and we found that it gives the same results as the more robust fourth and fifth order Runge-Kutta solver (ode45), while being considerably faster. The simulation is run multiple times and the linear phase of the average q apical curve is used to calculate the bTJ permeability coefficient: where w model is the width of the model system (m), which replaces the area in this 2D model. Next, a first degree polynomial is fitted to the linear phase of the average q apical to obtain the slope. Since the relative area of the junctions in the epithelium is considered in Eq 4, P bTJ is already scaled for the bTJs in the whole epithelium.
The tTJ central tubes are modeled as static pores, and an equation similar to Eq 4 is used: where � tTJ is the relative area of the tricellular pores in the epithelium, H p (r m /r tTJ ) is pore hindrance factor, r tTJ is the tricellular pore radius (m), and h tTJ is the tricellular pore height (m) [69]. The relative area of the pores is calculated as where ρ tTJ is the density of the tricellular junctions in an epithelium (m -2 ). The equation for H p was also derived by Dechadilok & Deen [70] by fitting a polynomial to computational results as where λ = r m /r tTJ . This equation is more accurate than the much-used Renkin equation when λ is close to unity [70]. The total epithelial TJ permeability is finally calculated based on the parallel connection between the two pathways as The required properties of the permeating molecule for this model are the molecular radius and the aqueous diffusion coefficient. Since we use a PEG oligomer with a mass of 547 Da, an equation relating the mass of a PEG oligomer to its size is used [71]: where r m is in Å (0.1 nm) and M m is the molecular mass (Da). The aqueous diffusion coefficient is calculated with an empirical relationship derived by Avdeef [61]: The default simulation time for the bTJ model is two hours and the stochastic behavior is further averaged by running the simulations 512 times. We found that the results were not affected by further averaging. The simulations were run using the Finnish IT Center for Science's (CSC) Taito supercluster using parallel computing (nodes with two 12-core Intel Haswell E5-2690v3 processors running at 2.6 GHz and 128 GB of DRR4 memory operating at 2133 MHz). A Matlab implementation of the permeability model is given in S1 File.

TER model
In the TER model, the bTJ and tTJ components are again calculated separately and connected in parallel in the end. The bTJ resistance model is constructed as a network of dynamic resistors ( Fig 2B). Since we are interested in the resistance, the strand capacitance is ignored. For each current loop i (Fig 2B), the equation is of the form where R ij (t) is the time-dependent resistance of the section of strand that is shared by current loops i and j (O), I i is the current in loop i (A), and V s is the voltage applied by the external source in the outer loop (V). The strand dynamics are incorporated into the resistances R ij (t), making them analogous to the rate constants in the bTJ permeability model. Because ions can also pass through the intact strands, R ij (t) depends on both the intact strand and break resistances: where l ij , l break , and r ij (t) have been described in the bTJ permeability model, R strand is the intact strand resistance per strand length (Om), and R break is the resistance of a break (O). The break resistance is calculated with the equation where ρ em is the resistivity of the extracellular medium (Om) and A break is the area of the break in the strand (m 2 ), calculated as 2w TJ l break . Although the TER measurement is basically instantaneous, the bTJ resistance model is simulated for a long time to average the results. The current flowing in the outer loop I outer (A) is used to calculate the bTJ resistance at each time point with Ohm's law: where the factor w model /l cb scales the results for the whole epithelium. To solve the bTJ model, the linear system defined by Eq 15 is transformed to matrix form and solved using Matlab.
The pores in the tTJ tubes are modeled as a static and their resistance is calculated as For each simulation time point of the bTJ resistance model, the total TER is calculated based on the parallel connection as Finally, to obtain the average TER for the simulation, the time average is taken from the results. The simulation time is 10 6 seconds and the simulations were run using the CSC's Taito supercluster with serial computing. A Matlab implementation of the TER model is given in S2 File.

Parameter values
Here we describe the default values of the model parameters. The TJ structure in our model is described by bTJ compartment dimensions, strand number (n strand ), strand height (h strand ), TJ half-width (w TJ ), tricellular pore radius (r tTJ ), and tricellular pore height (h tTJ ). The TJ dynamics are described by break forming and sealing probabilities (p break and p seal , respectively), and break size (l break ).
Although there is a lot a variety in the TJ strand morphology [5,[72][73][74][75], we used one strand morphology since the main focus is the break dynamics. The chosen bTJ compartment width (w comp ) and height (h comp ) were both 100 nm. These values are in the range of the bTJ compartment sizes found in the freeze-fracture replicas [5,[72][73][74][75] as well as by Kaufmann et al. using super-resolution microscopy [76]. The horizontal number of the compartments in the simulated systems was set to 50 and the heights of the apical and basal compartments in the molecular permeability model were both set to 200 nm. Based on the strand numbers in MDCK monolayers (3-5 strands) [52,53,74,77], Caco-2 monolayers (4-5 strands) [78], and retinal pigment epithelium (4 strands) [39], the default strand number was set to n strand = 4.
The value of h strand = 6 nm was based on the electron microscopy of TJ freeze-fracture replicas and the TJ strand architecture model by Suzuki et al. [14]. TJ half-width was chosen as w TJ = 4 nm, estimated based on the architecture model by Suzuki et al. [14] and transmission electron microscope images [54,79,80]. The dimensions of the tricellular pores, r tTJ = 5 nm and h tTJ = 1 μm, were taken from the measured values from freeze-fracture replicas [7,16].
The dynamic parameters of the strand dynamics were more uncertain due to the lack of experimental data. The strand breaks were assumed to remain open on average around 30 seconds based on the time scale of the dynamics in the transfected fibroblasts [46,50]. This led to the break sealing probability of p seal = 0.033 s -1 . The break size was approximated based on the figures and videos by Sasaki et al. [46] and by Van Itallie et al. [50], leading to a value of l break = 20 nm. This value was also used to quantify breaks by Rosenthal and coworkers [52,53]. The break forming probability (p break ) is fitted in the Results section based on the literature data. Due to the time scale of the dynamics and computational limitations, we restricted the possible state changes in the strands to occur every second.
The basal compartment concentration and the voltage of the external source used to measure the resistance are scaling parameters and do not affect the results. The chosen values were c basal = 1 mM and V s = 1 V, respectively. Also, the resistive properties of the breaks and the strands are needed. The resistivity of the extracellular medium required for the breaks was ρ em = 0.537 Om [16]. The value of the strand resistance (R strand ) depends on the epithelium, and is fitted in the Results section.
The cell boundary length per epithelial area (l cb ) and tTJ pore density (ρ tTJ ) are also highly dependent on the epithelium. They were determined using ImageJ Fiji [81,82] from the immunofluorescence microscopy images illustrating the cell-cell junctions in the studies our models were fitted for, and the values are described in the Results section. The only unknown parameters values were p break (both permeability and TER models) and R strand (TER model). These values are found by iteratively fitting the models to the experimental data. The model parameters described here are summarized in Table 1.

Model fitting and the origin of the leak pathway
The models were used to study the roles of tricellular junctions and bicellular strand dynamics in the leak pathway for the epithelial molecular permeability and TER. The models were fitted to the experimental data by varying the values of the break forming probability (p break ) and the TJ strand resistance (R strand ). Since the cell boundary length (l cb ) and tTJ density (ρ tTJ ) had a strong impact on the simulation results (see Parameter sensitivity analysis), we only used experimental data that included immunofluorescence microscopy images showing the cell-cell junctions. Therefore, unfortunately, the PEG profiling studies by Watson et al. [31,32], Van Itallie et al. [35], and Linnankoski et al. [36] had to be excluded from our model fitting.
First, the permeability model was fitted to the experimental data of 547-Da PEG oligomer permeation, since this molecule utilizes the leak pathway and it was used in the suitable studies [33,34,54,83]. The fitting was done by iteratively changing the value of p break (with the accuracy of 0.001 μm -1 s -1 ) and comparing the simulation result with the experimental result. The value of p break for MDCK C7 was calculated rather than fitted since the tTJ pores were enough to form the leak pathway for this epithelium, and thus making the fitting impossible. The value was iteratively calculated with Eq 3 using the chosen break sealing probability and the average amount of breaks per strand length for high-TER MDCK [52,53]. Next, the TER model was fitted using the obtained values of p break to iteratively find the values of R strand (with the accuracy of 0.01 GO μm). The experimental data used to fit the model, the values of experimental permeability and TER, the values of l cb and ρ tTJ , as well as the fitting results are shown in Table 2.

Structural dynamics in tight junctions
The cell sizes in MDCK II monolayers, as indicated by the cell boundary length per area (l cb ) and and the tricellular pore density (ρ tTJ ) in Table 2, varied greatly between the measurements. However, they were on average the largest of the fitted epithelia. The cells in Caco-2 monolayer were the smallest and in MDCK C7 monolayer between the two extremes. Surprisingly, the experimental permeability of high-resistance MDCK C7 was higher than that of the low-resistance MDCK IIb, which is most likely explained by the difference in cell size.
Although there was some variation in the values of p break for MDCK II, they are similar to each other having a mean (±SD) of 0.036 (±0.007). This is especially interesting when considering the great variability in the cell size. As for the TER, the variation was higher, with a mean (±SD) of 0.36 (±0.12) GO μm. The values of both p break and R strand were found to differ significantly for MDCK C7. Its values of p break and R strand were 7.2 times lower and over 30 times larger, respectively, compared with those of the average MDCK II. The properties of Caco-2 were a combination of the two MDCK strains: While the value of p break was similar to that of the MDCK II, R strand was closer to the value of MDCK C7. The resistance of a single 20-nm break in the strands was R break = 0.2 GO. The resistances of the same length of strand for average MDCK II, MDCK C7, and Caco-2 were 18 GO, 531 GO, and 383 GO, respectively. Thus, the strands had 90 to 2670 times higher resistance compared with the breaks.
To further check the validity of our results, we used the MDCK C7 p break and R strand values to simulate the TER for the original MDCK C7 epithelia (measured TER 5650 O cm 2 ) [84] by changing the cell size. Unlike with the other results considered here, the figures in [84] did not allow rigorous determination of the cell size properties, and therefore the cell radius was estimated to be between 15 and 20 μm (Fig 1A in [84]). Assuming perfect hexagonal cell array, we calculated the values of l cb and ρ tTJ to be between 0.050-0.067 μm -1 and 0.001-0.003 μm -2 , respectively. The obtained TER values for these cell radii of 15 and 20 μm were 4820 and 7310 O cm 2 , respectively. This indicated that the difference in cell size explains the difference in TER between the two experimental MDCK C7 results. Fig 3 shows the relative contributions of the two assumed leak pathway components-the tTJ pores and the bTJ strand dynamics-to the total leak pathway. MDCK C7 was the only epithelia considered here whose permeability was dominated by the tTJ pathway. Since the leak pathway was originally defined by the permeability, it can be said that the MDCK C7 leak pathway was fully formed by the tTJ pores. The role of bTJ dynamics was more prevalent, but variable, for the MDCK II permeability, having a mean (±SD) of 60.1 (±21.7)%. The bTJ dynamics was also the main pathway for the Caco-2 permeability. As for the TER, the role of the tTJ pores was insignificant for the MDCK II with a mean (±SD) relative role of 2.0 (±1.7)% between the four measurements. In both MDCK C7 and Caco-2, the tTJ pathway formed approximately half of the resistance of the epithelium.

Simulating experimentally-induced changes in the TJs
Next, we investigated how the developed models can recapitulate disturbances or changes in the TJ proteins, based on the studies of ZO-1 knockdown in MDCK II by Van Itallie et al. [34] and double ZO-1/2 knockdown in MDCK II by Fanning et al. [54]. Since both suggested that the observed changes in the barrier properties caused by these knockdowns were a result of decreased strand stability, we fitted our models to these results by varying the break forming probability p break . The TER model also required fitting of R strand . The results of the fitting are shown in Table 2  According to our results, the knockdown of ZO-1 led to a 62% increase in p break and to a 41% increase in R strand . In addition, the decreased bTJ tightness resulted in a 65% increase in the relative role of the bTJ pathway. The effect of the ZO-1/2 double knockdown was larger; it caused a 121% increase in p break and a 70% increase in R strand . The increase in the relative role of bTJ caused by the double knockdown (44%) was smaller than that of the ZO-1 knockdown due to the higher original contribution of bTJ in MDCK IIc. The changes in the relative roles of the pathways in TER were insignificant for both the single and double knockdown.

The effect of strand number on permeability and TER
Next, we studied the effect of the number of strands on the barrier properties by changing the strand number (n strand ) for the average MDCK II and MDCK C7 epithelia. These monolayers were chosen to illustrate the effect of n strand for different levels of strand dynamics and strand Structural dynamics in tight junctions resistances. Although these epithelia do not necessarily manifest varying numbers of strands, they provide two systems with different properties to base our simulations on. We simulated the model with n strand = [2, 6] for both permeability and TER, and the results for these simulations are shown in Fig 4A and 4B, respectively. To remove the impact of the cell size from the comparison, the simulations were run with the mean values of cell boundary length per area (l cb = 0.282 μm -1 ) and tricellular TJ pore density (ρ tTJ = 0.049 μm -2 ) of all the MDCK data included here (2 and C7).
Naturally, the permeability decreased as n strand increased (Fig 4A). In addition, the increase in n strand also led to an increase in the relative role of the tTJ pathway of the total permeability. This growing importance of the tTJ pathway led to saturation of the permeability at approximately 6 and 3 strands for MDCK II and C7, respectively. In contrast, TER grew when n strand increased ( Fig 4B). However, similarly to permeability, the significance of the tTJ pores increased with n strand . TER also showed the saturation behavior, but the saturation occurred past the simulated strand numbers for both of the MDCK strains. Moreover, the scale of the changes caused by the varying n strand were considerable larger for MDCK II in both permeability and TER. Also, the largest difference in permeability relative to the 4-strand standard system was almost two orders of magnitude compared with the under one order of magnitude for TER.
The raw, unaveraged simulation data indicating the time course of the simulations (Fig 4C  and 4D) showed the different behavior in the MDCK II and C7 monolayers for both permeability and TER. Full opening events, in which there was an open pathway through the strand network via the breaks, are indicated in the TER results by the sharp downward spikes. The spikes disappeared altogether at 6 strands for MDCK II and at 3 strands for MDCK C7. These events were not always clear in the permeability results, since the sharp steps in q apical may be a result of molecules stored into the small compartments released into the apical compartment. For example, as indicated by the TER results of 2-strand MDCK II, there was at least one full opening present around half of the time. However, no sharp steps were seen in q apical in any of the simulations shown. In contrast, there were multiple, minuscule scale steps in q apical e.g. for 6-strand MDCK II.
As n strand increased, the slopes for the bulk of the permeability simulations decreased for both MDCK II and C7 (Fig 4C and 4D). However, the increase in n strand in MDCK II led to more variable simulation results as well as to an increase in the number of the visibly different q apical curves with full openings. Interestingly, the behavior of MDCK C7 was different; generally, when n strand increased, the variability in the results decreased. This was most likely due to the lack of full opening events. Furthermore, the clearly distinct q apical curves in the 3-strand MDCK C7 differed considerably more from the bulk of the simulations compared with other simulated systems.
Thus, the simulation raw data showed a biphasic behavior for bTJ permeability as a function of TJ tightness, defined by both the strand number and the level of strand dynamics. With low strand numbers and high values of break forming probability, the bulk of permeability through the strand network passed via the full opening events. This resulted in low variance between the individual simulations, as shown, e.g., by the 2-strand MDCK II (Fig 4C). With high strand numbers and low values of levels of strand dynamics, however, the bulk of the permeability occurred via the step-by-step diffusion through the compartment network. Again, this resulted in low variance between the individual simulations, as shown, e.g., by the 5-and 6-strand MDCK C7 (Fig 4C). Between these two extremes, there was a transition zone where the variance between the simulations was higher.
The permeability model was simulated with the equilibrium state of the system as the initial condition. These states for each of the compartment rows were found by running the models with zero initial conditions for a long time. The relative equilibrium concentrations compared to the basal compartment are shown in Fig 5 for systems with 2-6 strands. The parameters affecting the rate constants or the magnitude of p break had no effect on these values; they only defined how fast the linear phase was reached. Interestingly, while the differences between the small TJ compartments were approximately linear, the equilibrium concentration change between the basal compartment and the bottom small compartment row as well as the top small compartment row and the apical compartment showed nonlinear discontinuities.
In addition, although the lag phase was not included in the simulations, we calculated the length of this phase for each of the bTJ permeability simulations to describe the time it takes for the permeating molecules to pass through the TJs. This was done by running the simulation with zero initial values for the small bTJ compartments and by extrapolating the linear phase of q apical backwards to determine its intersection with the time axis. Longer than normal simulation times were used in some cases to obtain a linear phase of sufficient length. The lag times for both of the MDCK strains are shown in Table 3. As expected, the lag time grew as the strand number increased and as the break forming probability decreased. Interestingly, the lag  time for 5-strand MDCK II and 2-strand MDCK C7 were close to each other, although having a large difference in the actual permeability coefficients. The same biphasic behavior can be seen in the lag times when the barrier became tighter. The lag times grew slowly for MDCK II as n strand was increased. However, with the higher strand numbers for MDCK II and for all the results for MDCK C7, the increase in n strand led to considerably larger changes in the lag times.

Comparison with steady-state models
To test if the results produced by the dynamic bTJ models presented here could be reached with simpler methods, we created steady-state (SS) bTJ models that assumed a static system with an average number of breaks per strand for both barrier properties. In the SS bTJ permeability model, the compartment rows were reduced into a single large compartment between the strands, since the compartments in the same row could be assumed to be in equilibrium.
The SS bTJ resistance model similarly assumed to only contain horizontal strands, and thus simplifying the model to a series connection of identical strands. The number of static open breaks for both SS models was defined from Eq 3. For comparison, we ran the simulations for the varying strand number for both MDCK II and C7 presented in the previous section, and the results of the comparison are shown in Fig 6. It can be clearly seen that the SS models were not able to produce comparable results for bTJ permeability nor for bTJ resistance. The permeabilities predicted by the SS model were well above those produced by our dynamic model (Fig 6A). Moreover, the SS permeabilities showed very little change overall when the strand number was increased from 2 to 6. As for the bTJ resistance, the values produced by the SS model were approximately one and two orders of magnitude lower, respectively, than those from our dynamic model for MDCK II and C7 (Fig 6B).

Parameter sensitivity analysis
Finally, we performed a sensitivity analysis on certain model parameters by individually altering their values ±25%. The chosen parameters for both models were break forming and sealing probabilities (p break and p seal , respectively), break size (l break ), cell boundary length per area (l cb ), and tTJ pore density (ρ tTJ ). Strand resistance (R strand ) was additionally included for the Structural dynamics in tight junctions TER model. Parameters l cb and ρ tTJ both depend on the cell size, and thus are not independent from each other. However, by changing them individually we could observe the relative roles of these parameters. We ran the analysis for both the average MDCK II and the MDCK C7 using the same default values of l cb (0.282 μm -1 ) and ρ tTJ (0.049 μm -2 ) as with the strand number simulations. The results comparing the sensitivity simulations with the standard simulations are shown in Fig 7. The permeability of MDCK II was very sensitive to the changes in the strand dynamics, since alterations in the break probabilities led to large changes in the permeability. The variance of these parameters had a significantly smaller effect on the MDCK II TER. While the alterations in l break affected the MDCK II permeability, they had no effect on TER. Alterations in the values of l cb and ρ tTJ had similar levels of influence on the MDCK II permeability. The MDCK II TER was unchanged by the variance of ρ tTJ . However, it was affected more by the variance of l cb than the permeability. Finally, the changes in TER were equal to the alterations in R strand for MDCK II, indicating direct proportionality.
The results were extremely different for MDCK C7. Due to the lower level of strand dynamics, the alterations in the parameters describing the breaks (p break , p seal , and l break ) had no effect on either permeability or TER. This was also the case with l cb for permeability. On the other hand, the results indicate direct dependence of permeability on ρ tTJ for MDCK C7. The impacts of both l cb and ρ tTJ were similar for TER, but smaller than the parameter value variations. Finally, the influence of R strand for MDCK C7 was smaller than for MDCK II. Results of the parameter sensitivity analysis. We varied the values of the chosen parameters by ±25% and both the permeability and TER simulation results were compared with the normal system. The analysis was conducted for average MDCK II permeability (A) and TER (C) as well as MDCK C7 permeability (B) and TER (D). p break , break forming probability; p seal , break sealing probability; l break , break size; l cb , strand length per area; ρ tTJ , tricellular TJ pore density; R strand , strand resistance. https://doi.org/10.1371/journal.pone.0214876.g007

Discussion
Tight junctions (TJs) are an indispensable part of the epithelia that form the barriers between many of the body's compartments, and yet not enough is known about their structure or structural dynamics. In this work, we have developed a computational model of the dynamic TJ structure to study the origin and the properties of the leak pathway-the nonspecific permeation pathway through the TJs. This was done by simulating the epithelial molecular permeability of a PEG oligomer and transepithelial electrical resistance (TER) with the same structural strand dynamics for different epithelial monolayers and scenarios. The model combines the current knowledge and theory of the dynamic TJ structure into a computational framework.
There are experimental findings that attest to both candidates for the leak pathway: the tricellular junction pores and the bicellular strand dynamics. Krug et al. [16] observed that 3-kDa dextran mainly diffuses through the tTJ pores, and through the bTJs to a lesser extent. While our model did not extend to macromolecule permeation, the limited size of the strand breaks would result in the tTJ pores being the main pathway for the molecules of this size. Krug et al. also calculated the role of the tTJ pores to be minuscule for ion permeation due to their rarity [16]. At the moment, there is no direct evidence of the dynamic bTJ strand breaking and sealing events forming the leak pathway, but it has been theorized [20,22,33,47,48]. Although dynamic strand breaks have been observed in transfected fibroblasts [46,50] and static breaks in freeze-fracture images [16,[51][52][53], no strand breaks were detected by Weber et al. in their bTJ patch clamp measurements [56]. They discussed multiple reasons for this, including that the breaks may not be distinguishable from the patch seal break [56]. It is, however, also possible that the strong electrode seal between the pipette and the junctional membrane might stabilize the strands mechanically and therefore prevent the strand level dynamics.
While one of the assumptions in our models was that the leak pathway is formed by both the bTJ strand dynamics and tTJ pores, interestingly, the permeability leak pathway of the MDCK C7 was formed solely by the tTJ pores. Whereas the relative roles of two pathways in permeability varied greatly between our four MDCK II fitting results, they all showed that the tTJ pores by themselves were incapable of producing the measured permeability values. A possible cause for the variation between the MDCK II results is the cell culture times in the experiments. The reported respective culturing periods for the MDCK IIa, IIb, and IId were from 4 to 8 days, from 7 to 10 days, and 10 days [33,34,54]. There was no culture time directly reported for MDCK IIc [83]. This indicates that longer culturing periods might have led to a higher significance of the bTJ pathway and to a higher level of strand dynamics. However, since the amount of the data is limited, this might be a coincidence. Similar to MDCK II, approximately half of the permeability leak pathway was formed by the bTJ dynamics for Caco-2. As for the TER, the MDCK II was dominated by the bTJ pathway, as was also calculated by Krug et al. [16]. On the contrary, the impact of tTJ pores on the Caco-2 and the MDCK C7 was significant, due to the higher strand resistance.
Thus, as theorized by Liang & Weber [20], our model suggests that the tTJ pores and the bTJ strand dynamics may both contribute to the leak pathway with varying degrees. Moreover, the significance of these two alternatives was different for permeability and TER. The tTJ pores were the prominent permeation pathway only for epithelia with more stable strands, since the extremely rare breaks did not enable fast step-by-step permeation. The strand dynamics had a lesser role in determining the bTJ resistance, as the strand resistance, and thus their molecular composition was the dominating factor. The tTJ pores became more important only with the higher strand resistances. This two-component leak pathway is further supported by how changes in the expression of occludin and tricellulin, effect the macromolecular permeation.
Although it is not completely understood how, the expression of levels of occludin, mainly located in the bTJ strands, have been shown to regulate the leak pathway permeability [45,83,85,86]. On the other hand, Krug et al. [16,87] have shown that increased expression of tricellulin leads to a decreased permeability of >4-kDa dextrans and vice versa.
Although there was diversity in the cell sizes, the main differences between the epithelial barriers rose from the distinct levels of the bTJ strand dynamics and strand resistances. Of the epithelia considered here, Caco-2 was found to have the most dynamic strands. The obtained break forming probabilities for the MDCK II results were quite similar to each other, while MDCK C7 had extremely limited strand dynamics. Surprisingly, while the strands of Caco-2 were the most dynamic, its strand resistance was comparable to that of the MDCK C7. In addition, the large difference in measured TER between MDCK II and C7 was also observed in the obtained strand resistances, which is in line with the current understanding that the TER is mainly defined by the conductive properties of the claudins, especially claudin-2, found in the strands [33,74,88]. The strand resistances of every epithelia were significantly higher than the strand break resistance. However, due to the rarity of the breaks, especially full openings, the strand resistances largely defined the overall TER, as shown by the sensitivity analysis. The MDCK C7 TER measured by Van Itallie et al. [33] was considerably smaller than the originally measured values (460 vs. 5650 O cm 2 ) [84]. Our simulations show that this difference is explained by the cell size, since the C7 cells of Van Itallie et al. [33] are distinctly smaller than those of Gekle et al. [84], which leads to smaller length of bTJs and number of tTJ per area and thus higher TER. Further, although the claudin-2 dynamics model by Weber et al. [66] used a 3-strand model instead of the 4-strand used here and their description of the leak pathway differed from ours, there are surprising similarities. The calculated steady-state strand resistance defined for their claudin-2-transfected high-resistance MDCK I model equals to approximately 0.32 GO μm, which is close to our claudin-2-containing MDCK II values.
Our models showed what kind of changes in the dynamic structure could lead to the observed experimental changes in the TJ barrier properties. Van Itallie et al. [34] found that ZO-1 stabilizes the TJ barrier, and therefore we described the decrease in stability caused by the ZO-1 knockdown by an increase in the strand break forming probability. The permeability results of the double ZO-1/2 knockdown study by Fanning et al. [54] were similarly replicated by changing this parameter. The higher change in the break forming probability in the double knockdown is in line with the observations that the two ZOs have redundant roles [54,89], and thus knocking out both of them should decrease the strand stability even more. The decreased stability caused by the lack of binding between claudin-2 and ZO-1 was recently visualized in fibroblasts by Van Itallie et al. [50]. Interestingly, they found no difference in the number of breaks with or without this binding. However, this result might have been caused by the nonepithelial model system. On the other hand, we did not consider possible changes in the bTJ strand number or morphology due to the knockdowns. Umeda et al. [89] showed that ZO-1 knockout/ZO-2 knockdown completely eliminated the bTJ strands, accompanied by extremely low TER compared with the control. However, we expected that the strand number was not greatly affected by the ZO-1 and double ZO-1/2 knockdowns, since only a minor change in TER was reported [34,54]. In addition, since we also had to increase the strand resistances to fit the TER model to these data, our results indicate that these knockdowns might also affect the strand resistance in a presently unknown manner.
The understanding regarding the role of the bTJ strand number in the properties of the epithelial barrier has changed over time. TER was originally found to depend exponentially on bTJ strand number [21]. This was hypothesized to arise from transient pores in the compartmentalized bTJ strand network [21,90]. However, this idea has been refuted since TER is now known to mainly depend on the TJ claudin composition [5,74,88]. Our results agree that the dependence of TER on the strand number is not straightforward and show a complex dependence on strand resistance that is defined by the protein composition, the strand number, as well as the level of structural dynamics. The saturation towards higher strand numbers comes from the increasing relative role of the static tTJ pathway. Moreover, the results of TER as a function of time show the immensely varying behavior of the resistance in the epithelia with different levels of strand dynamics and numbers.
To the best of our knowledge, there have been no experimental studies that directly investigate the effect of strand number on molecular permeability. Colegio et al. [75] showed that an increase in strand number caused by the increased claudin-2 and -4 expression had no effect on the permeability of mannitol that diffuses via the leak pathway. However, based on their freeze-fracture images [75], the strand numbers were in the range of the permeability saturation shown in our simulations, and thus mannitol permeability should remain unchanged, indicating agreement with our results. The saturation in our results was caused by the tTJ pores, as the permeability of the bTJ pathway become lower than that of the tTJ pathway. The unaveraged permeability simulation raw data revealed vastly different permeability behavior depending on the strand numbers and levels of breaking dynamics. The observed biphasic behavior in the permeability simulations is most likely a property of the dynamic network system. In the high-permeability side of this behavior, the change in the number of the full opening events in the strand network caused by alterations in the strand number or level of dynamics led to large changes in permeability. However, since the low-permeability end of this behavior depends on the step-by-step diffusion between the compartments, the alterations in the strand number or dynamics have a lesser effect on the permeability.
Overall, our results concerning the effect of strand number, the ZO knockdowns, as well as the sensitivity analysis showed that changes in the bTJ structure or dynamics lead to considerably larger changes in the molecular permeability compared with the TER. Thus, measuring only the TER might hide important unseen changes in molecular permeability and in the barrier properties. The TER is a good indicator of epithelial condition and is straightforward to measure, but it contains a lot of uncertainty due to the differences in measurement setup and conditions. Our results strengthen the idea that both the molecular permeability and the TER are needed to properly define the TJ as a barrier [91,92].
Our comparison with the steady-state models showed that they are not able to reproduce the same behavior as our dynamic models. Since there were constant breaks in the steady-state model, the changes in the strand number or the number of breaks had only a minor effect on the permeability and the TER. Nevertheless, our models had their limitations concerning both the chosen parameter values and the geometry. Because of the lack of rigorous experimental data, we were forced to estimate and fit the probabilities of the strand dynamics and build the model partially with assumptions and hypotheses available in the literature. Moreover, the TJ strand morphology is very diverse and heterogeneous, the number of strands varies within an epithelium [39,55,77,78,93], and the strands typically become tighter towards the apical direction [13]. Also, the expression of tricellulin differs between the epithelia, but the effects of its expression level on the TER and the permeability of molecules with similar size to 547-Da PEG were minor [16]. Further, many molecules permeate the epithelia through the cells using active transport processes, which are difficult to model due to their specificity. However, the PEG oligomers are hydrophilic [40], indicating that they mainly diffuse via the TJs. Thus, our model reflects the PEG-based permeability measurements. In addition, our model lacks the claudin-2 pore dynamics that have previously been observed [56,66]. However, we opted to exclude these dynamics, since their effect has not been characterized for noncharged molecules. Finally, although there is evidence that the TJs form the main conduction pathway [27,28], Günzel et al. [26] calculated that for many epithelia the transcellular resistance affects the overall TER or is even lower than the paracellular resistance. This indicates that our assumption that TJs define the TER will not work for every epithelium and condition.
We demonstrated that our dynamic structure models can reproduce the basic TJ barrier behavior. With the simplified structure and dynamics, the models enable the comparison between the molecular permeability and the TER under the same dynamic context at the level of the TJs. We showed that the TJ strand breaking dynamics can drastically alter both of these barrier properties independently of each other, highlighting the importance of both measures for characterizing the epithelial barrier. Furthermore, our results indicated that the leak pathway may be formed both by the tTJ pores and the bTJ strand dynamics with varying degrees, but differently for the permeability and the TER. Our models create a good methodological framework that can be used to integrate knowledge on TJ structure, parametrize experimental measurements, and produce hypotheses that can be studied experimentally. Refined versions of the models could include a more realistic strand network as well as inhomogeneities in the strands. This would provide tools to study how diseases that affect the TJ structure alter the properties of the epithelial barrier.