Mechanism of interdigitation formation at apical boundary of MDCK cell

Summary It has been reported that the MDCK cell tight junction shows stochastic fluctuation and forms the interdigitation structure, but the mechanism of the pattern formation remains to be elucidated. In the present study, we first quantified the shape of the cell-cell boundary at the initial phase of pattern formation. We found that the Fourier transform of the boundary shape shows linearity in the log-log plot, indicating the existence of scaling. Next, we tested several working hypotheses and found that the Edwards-Wilkinson equation, which consists of stochastic movement and boundary shortening, can reproduce the scaling property. Next, we examined the molecular nature of stochastic movement and found that myosin light chain puncta may be responsible. Quantification of boundary shortening indicates that mechanical property change may also play some role. Physiological meaning and scaling properties of the cell-cell boundary are discussed.


MDCK cell tight junction shows a ragged interface with temporal fluctuations
Quantification of the ragged cell-cell interface revealed scaling of the boundary shape Stochastic movement and boundary shortening can explain the observed scaling Myosin puncta may be responsible for the stochastic movement of cell-cell boundaries

INTRODUCTION
Various cells show interdigitated cell boundaries. For example, plant leaf epidermal cells show a winding pattern of cell walls, 1,2 which reduces the mechanical stress exerted on the cell wall. 3 Kidney podocytes show intricated interdigitation pattern, 4 which works as a filter to generate urine. 5 Endothelial cells show interdigitated cell boundaries during collective cell migration. 6 High glucose stimulus induces interdigitation of tight junction and barrier dysfunction of small intestine epithelium. 7 Various other transporting epithelia, including renal tubules, 8 avian salt glands, 9 elasmobranch rectal glands, 10,11 gill epithelia, [12][13][14] and yolk-sac membrane in fishes, 15 have been reported to show interdigitated cell junctions. All of these are ''leaky'' epithelia with high paracellular conductance, and the interdigitated pattern has been proposed to enhance paracellular transport. However, the mechanisms underlying the formation of the interdigitated cell boundary and its physiological significance are not fully understood.
Several mechanisms have been proposed to explain the interdigitation pattern formation at cell-cell junctions. In plant leaf epidermal cells, cell wall synthesis and degradation by ROP can be a primary mechanism for interdigitation. 16 We used this mechanism and established a theoretical model 1 based on our previous model of skull suture interdigitation formation 17,18 to reproduce this interdigitation formation. In addition, we pointed out the possibility that the buckling instability may also play a role in the pattern formation. 2 In the case of endothelial cells, the filopodia are inserted into the cytoplasm of the preceding cells, which appears as cell-cell interdigitation. 6 Madin-Darby canine kidney (MDCK) cells are a model mammalian cell line that is utilized for epithelial pattern formation study 19 and assaying epithelial barrier. 20 It is reported that the interdigitated pattern is generated in the tight junctions of MDCK cell sheets ( Figures 1A and 1B). 21 The interdigitation is limited to the apical tight junction region, so the structure cannot be observed with normal brightfield observations. The junction is not a static structure but shows fluctuation during the pattern formation. 21 Kidney epithelium at the loop of Henle also shows extensive interdigitation in vivo ( Figure S6B), 8 suggesting this interdigitation may reflect the morphology in vivo. However, the pattern formation mechanism of this apical interdigitation remains to be elucidated.
In the present study, we first quantified the cell interdigitation dynamics in the MDCK cell line. Next, we quantitatively characterized the pattern and found that the Fourier transform of the interdigitation pattern showed linearity in the log-log plot, indicating the scaling property of the system. Then we chose four working hypotheses for spontaneous pattern formation and verified the models by the quantitative measurements. Among these hypotheses, we chose the Edward-Wilkinson model which consists of stochastic fluctuation and shortening of the boundary that can generate the observed scaling property. Then we searched for possible mechanisms of stochastic fluctuation by immunohistochemistry, live imaging, and inhibitor assays and showed that myosin puncta might be responsible for the random movement of cell-cell junctions. Finally, we separated and quantified the effects of boundary shortening and random fluctuation using the image processing technique.

Observation of interdigitation in MDCK cells
First, we used ZO-1-EGFP MDCK cells to observe the interdigitation of cell-cell junctions. We could observe the formation of interdigitation in MDCK cells ( Figure 1A). To observe the tight junction morphology at the ultrastructural level, we observed the freeze-fracture replica of the apical plasma membrane of MDCK II cells. The interdigitation of the tight junctions was visible in the replica that spanned the apical membrane, consistent with the immunofluorescent staining ( Figure 1B). This interdigitation pattern was also observed in kidney tubules in vivo. Immunohistochemistry of tight junction and thick ascending limb markers showed that NKCC2-positive tubules have strongly curved tight junctions ( Figure S6B).

Live cell imaging of cell boundary dynamics
Next, we undertook long-term live imaging of MDCK cell boundaries by avoiding the dome (blister) formation. MDCK cells were cultivated under the standard protocol. Because the cells retain the characteristic to transfer iScience Article fluid, we frequently observed the dome formation and detachment of the epithelial sheet from the culture dish. 22 To minimize this effect, we cultivated the MDCK cell at the bottom of the culture insert and observed the dynamics from the bottom using a glass-bottom dish and a plastic spacer (Figures S1A-S1C). The cells grew and started showing interdigitation without forming domes, which enabled the detailed observation.
We observed the formation of interdigitation at the cell-cell junctions after one week of cultivation. Live imaging of the pattern formation process showed that the interdigitation structure shows stochastic movement ( Figures 1C-1F), as reported in previous studies. 21 Quantification of the interdigitation shape Next, we quantified the shape of the interdigitation. Cell boundaries of MDCK cells were obtained, and the boundary line shapes were converted to one-dimensional functions hðxÞ (Figures S1D-S1F). Then, the power spectrum ( Figure S1G), tortuosity ( Figure S1H), and the maximal amplitude ( Figure S1I) of each timepoint were obtained.
We found that Cj b hðkÞj 2 D showed linearity in the log-log plot, which indicated the existence of scaling (Figure 2C). We used 133 and 92 segments for day 0 and day 10 samples, respectively. Scaling exponents are 2.18 and 2.19 for day 0 and day 10, respectively. In addition, the linear distribution shifted upwards as the formation of interdigitation proceeded from day 0 to day 10 ( Figure 2C). This may imply that the complexity of the boundary shape increases (Figures 2A and 2B). We also obtained tortuosity and the maximum  iScience Article amplitude. Both of these indexes increased with time, indicating that the boundary became more complex ( Figures 2D and 2E).

Models of pattern formation that may explain the observed scaling
To understand the mechanism of pattern formation that generates scaling observed in Figure 2C, we chose four working hypotheses that may explain the characteristic scaling of the interdigitation pattern.

White noise model
One obvious candidate model of the pattern formation is the white noise, in which the membrane moves according to the random fluctuation ( Figure 3A). The governing equation is: v vt hðx; tÞ = hðx; tÞ: (Equation 1) hðx; tÞ represents white noise. We defined the spatial average of the noise term as 0 because the boundary of two equivalent cells should not have any preferential direction of movement. We defined the mean and variance of the noise term h as 0 and D, respectively (hðx; tÞ = 0, hðx; tÞ 2 = D). h represents the spatial average of h. In this model, all the membrane fragments move independently. As a result of the numerical simulation, we obtained a ragged boundary shape, which resembled the observed pattern ( Figure 3E).
However, the sample mean of Fourier transform of the pattern Cj b hðkÞjD showed a flat pattern ( Figure 3I, a solid line), which was different from the observed pattern ( Figure 3I, dashed line).

Edwards-Wilkinson model
Another possible candidate is the Edwards-Wilkinson model ( Figure 3B). 23  hðx; tÞ represents the stochastic term defined in the previous section, and d h is the strength of surface tension that minimizes the length of the boundary.
Numerical simulation of the model generated a randomly fluctuating boundary ( Figure 3F) that looked similar to the observed boundary shape (Figures 2A and 2B). Fourier transform of the pattern Cj b hðkÞjD showed linearity ( Figures 3J, a solid line), which was the same as the observed pattern ( Figure 3J, a dashed line). Moreover, the gradient of the line of the simulation pattern was identical to that of the experimental observation, whichfurthersupported the model. It is known that the gradient of the power spectrum is À 2 at the steady state of the Edwards-Wilkinson Model (supplemental information S2).

Turing model
One of the most famous working hypotheses of pattern formation is the Turing model ( Figure 3C) 24  In this case, the shape hðx; tÞ is determined by the distribution of some diffusing signaling molecules uðx; tÞ and/or vðx; tÞ. f ðu; vÞ and gðu; vÞ represent the interaction between u and v, and d u and d v represent the diffusion coefficients of u and v, respectively. The shape of the boundary is determined downstream of these molecules, for example, v vt hðx; tÞ = eðuðx; tÞ; hðx; tÞÞ.  iScience Article Numerical simulation of the model produced a regular interface shape, which looked different from the actual pattern ( Figure 3G). Fourier transformation of the generated pattern ( Figure 3K, a solid line) showed a sharp peak representing the characteristic length of the generated pattern. Overall, the pattern generated by the model was qualitatively different from the observed pattern (Figures 3K, a dashed line).

Buckling model
Another well-known mechanism of pattern formation is buckling instability ( Figure 3D), in which boundary length is relatively longer than the cell volume constriction, and has been implicated in the boundary interdigitation of Drosophila aminoserosa cells. 25 The linear part of the buckling model, which represents the onset of pattern formation, is described as follows (supplemental information S1, Figure S3) c is the elongation of the cell boundary, and d is the mechanical characteristics of the elastic bar, which tends to make the curved boundary straight.
Numerical simulation of the model also showed a relatively regular periodic pattern ( Figure 3H). In addition, the power spectrum of the simulation result ( Figure 3L, a solid line) was very different from the observed pattern ( Figure 3L, a dashed line).

Parameter change to reproduce upper shift of power spectrum
In the previous section, we showed that the Edwards-Wilkinson model is the most plausible working hypothesis that explains the generated pattern. Next, we predicted from the model which parameter change can reproduce the upper shift of j b hðk; tÞj 2 in the Edwards-Wilkinson model. Either decreasing the diffusion coefficient d h ( Figures 3M and 3N) or increasing the noise term hðx; tÞ (Figures 3M and 3O) can reproduce the upper shift of j b hðk; tÞj 2 ( Figure 3P). It is impossible to discriminate whether decreasing the diffusion coefficient (= softening of the boundary structure) or increasing the noise term (= increase of stochastic movement of the boundary) is responsible for the upper shift.

Mechanism of random fluctuation of the cell boundary Distribution of MHC-B and membrane tension
We observed the role of actomyosin in the boundary interdigitation. We used immunohistochemistry to compare the formation of interdigitation (ZO-1) and actomyosin system (MHC-B). We also used anti-vinculin and a18 antibodies to detect the tension applied to cell junctions. It is known that a-catenin changes its conformation by the application of tension ( Figure S2A). 26 On day 1, interdigitation was not obvious, and myosin distribution was not detectable. As the interdigitation developed on days 2 and 5, we observed the formation of myosin puncta ( Figure 4A, red). The images are captured using the same conditions, and the puncta are highlighted by using a constant threshold value. The number of puncta increased at day2 and day 5 compared to day 1 ( Figure 4B, ANOVA, p < 0:01). At the same time, vinculin staining at the cell boundaries became stronger, indicating the increased tension of the boundary. Triple staining of ZO-1, a18, and MHC ( Figure 4C) showed that the interdigitated boundary showed uneven tension distribution and MHC-B puncta. We undertook immunohistochemistry to further observe the myosin distribution ( Figure 4D). Complex myosin fibrous networks run through apical cytoplasm (white arrows), and these myosin fibers are focally attached to the junctions, as a myosin puncta (white arrowheads). We also observed the phosphorylation of myosin using ppT18S19, which stains fully active myosin ( Figure 4E). 27,28 We could also observe network-like structure (white arrows) and puncta-like staining at the junction of myosin fiber and ZO-1 (white arrowheads).

Role of actomyosin in interdigitation formation
To directly examine the mechanism of random fluctuation of the cell boundary, we observed the effect of various inhibitors of actomyosin activity on the interdigitation formation. iScience Article At first, we examined the role of myosin in interdigitation formation. Blebbistatin treatment, which inhibits the function of active myosin ( Figure S2B), causes the disappearance of the boundary interdigitation, and the cell boundary became completely straight. 21 We reproduced ( Figures 5A and 5B) and quantified the result ( Figures 5C-5E). The log-log plot of the amplitude of each wavenumber component showed a lower shift by the blebbistatin treatment ( Figure 5C). We used 718 and 687 segments for control and Blebisttatin samples, and the scaling exponents are À 2:18 and À 2:25 respectively. Tortuosity and amplitude of blebbistatin-treated samples were decreased ( Figures 5D and 5E). Statistically significant differences were detected (Mann-Whitney test, p < 0:01).
Y-27632, which blocks the effect of ROCK and inhibits the conversion of inactive myosin to active myosin ( Figure S2B), showed slightly different behavior ( Figures 5F-5J). The treated cells showed less interdigitation, but the boundary became a large arc, not a straight line ( Figure 5G, white arrows). As a result, we could observe a change of Cj b hðk; tÞj 2 D gradient (2.26 and 2.39) ( Figure 5H) and a reduction of tortuosity ( Figure 5I), but the maximal amplitude was not changed ( Figure 5J. n = 56 and 30 for control and Y-27632 samples respectively). Statistically significant change of tortuosity was detected between control and Y-27632treated groups (Mann-Whitney test, p < 0:01).
To examine the effect of myosin inhibitors on the distribution of phosphorylated myosin and myosin, we undertook immunohistochemistry after the myosin inhibitor treatment ( Figure 6). We used an antibody to detect myosin phosphorylated at Thr18 and Ser19, representing the fully active myosin. 27,28 The control sample shows prominent interdigitation of ZO-1 structure and complex network of myosin fiber ( Figure 6A, arrows), fully phosphorylated puncta ( Figure 6B, arrowheads) and strong uneven a18 spots ( Figure 6C, arrows). Blebbistatin treatment reduced the apical myosin network ( Figure 6D), but phosphorylation itself is not reduced ( Figure 6E). a18 expression was considerably reduced ( Figure 6F). Y-27632 treatment diminished the myosin network ( Figure 6G) and reduced the phosphorylated myosin light chain in the apical myosin network ( Figure 6H). Y-27632 also resulted in a relatively even a-18 distribution compared to the control ( Figure 6I). The effect of myosin inhibitors on the distribution of puncta and tension can be interpreted as follows: blebbistatin is an inhibitor of myosin II ATPase activity ( Figure S2B), so it does not interfere with the phosphorylation of myosin ( Figure 6E) but the tension is reduced by the relaxation of actomyosin microfilaments ( Figure 6F). Y-27632 is a Rho kinase inhibitor ( Figure S2B), which directly affects the phosphorylation status of the myosin ( Figure 6H). 29 We also examined the role of actin by inhibitors, but we could not obtain any statistically significant difference ( Figure S4). We tried SMIFH2, an inhibitor for Formin-mediated polymerization of unbranched actin, and CK666, which inhibits Arp2/3 complex that regulates polymerization of branched actin networks (Figure S2C). We quantified 777 control segments, 743 SMIFH segments and 785 CK666 segments. Scaling exponents are À 2:1 (Control), À 2:0 (SMIFH2) and À 2:1 (CK666). We could not observe a detectable difference ( Figures S4A-S4C), a shift of j b hðk; tÞj 2 ( Figure S4D), or a change of tortuosity ( Figure S4E) or maximal amplitude ( Figure S4F).

Dynamics of MRLC2 puncta
Next, we focused on nonmuscle myosin II, which is a motor protein that crosslinks actin filaments and generates contractile force. Nonmuscle myosin II was visualized by expression of myosin regulatory light chain 2 (MRLC2)-EGFP fusion protein. Firstly we cultivated 1:9 mixture of wildtype: MRLC2-EGFP cells. We could observe the boundary between these cells ( Figure 7A).
At first we confirmed the three-dimensional distribution of MRLC2 to select the focus plane to observe ( Figure 7A). We used a 1003 objective lens and the smallest pinhole size to avoid the fluorescence from different z planes. At the basal side of the cell, we observed strong fluorescence alongside the iScience Article putative actin stress fiber structure, which was not the target of our study. On the z plane of the apical side, we observed puncta structure within the cytoplasm and at the lateral membrane of the cells. We sometimes observed static, regularly distributed puncta at the apical junction level ( Figure 7A, a white circle). On the apical-most plane we observed many puncta at the apical membrane ( Figure 7A, white arrows).
Next, we observed the dynamics of the puncta. The puncta near the apical membrane showed dynamic movement in the kymograph (Figures 7B and 7C). Some of the puncta moved dynamically while others were static. Then we examined the spatial characteristics of the myosin distribution. Fourier transformation of the spatial and temporal pattern of the kymograph showed both spectra had cutoff frequencies ( Figures 7D and 7E, dashed lines). The amplitudes of waves with wavenumbers 1-4 were large, and those with higher wavenumbers were small. Therefore, the size of the smallest unit of the bright spot was Next, we observed whether the dynamic movement of puncta was correlated with the boundary movement by observing the boundary between wildtype and MRLC2-EGFP cells ( Figure 7F). We observed that bright iScience Article puncta sometimes appeared at the apicolateral membrane and the boundary movement could be concomitantly observed. The movement can be both forward ( Figures 7H, a solid arrow) and backward ( Figures 7G and 7H, dashed arrows). One round of push usually lasted around 10 min, and the membrane resumed its original position after the puncta disappeared ( Figures 7G and 7H). iScience Article

Estimation of d h of cell boundary
Next, we tried to separate the effect of d h from the experimental data to clarify whether d h was responsible for the upper shift of the line in Figure 2C. Detail of the method is described in supplemental information S4 and Figure S5. In short, if we assume that each point in a cell boundary obeys the Edwards-Wilkinson equation, the noise term h can be eliminated by taking an average of many datapoints.
We obtained d h at four timepoints (0, 30, 60, and 90 min) on day 0 and day 6, and compared the calculated diffusion coefficients (Figure 8). We detected a statistically significant decrease of d h between day 0 and day 6, indicating the change of mechanical characteristics also influenced the shift of the line in Figure 2C (Mann-Whitney test, p < 0:01).

Observation of cell boundary dynamics in vivo
One obvious question is whether this stochastic pattern formation process also takes place in vivo (Figure S6B). To examine this possibility, we used ZO-1-EGFP transgenic mice to observe whether the stochastic movement was observed in in vivo situation. Because the tight junction structure became disrupted swiftly after the preparation of ex vivo sample, we observed the cut surface of the kidney specimen immediately after the preparation. Because the region of tubules which had extensive interdigitation was limited, we minimized the search time of the thick ascending limbs of Henle by looking at the cortex-medulla boundary ( Figure S6A).
Against all these efforts, currently we only observed only a minor movement of the cell-cell boundary ( Figures S6C-S6E). Most of the boundaries remained unchanged ( Figure S6C), but in some cases we observed fluctuation of cell-cell boundary (Figures S6D and S6E).

Physiological roles of interdigitation formation
The physiological meaning of this apical interdigitation is yet to be elucidated. One hypothesis is that interdigitation may influence paracellular transport. 30,31 Among the three types of junctional complexes, 32 tight junctions regulate paracellular transport, and the claudin proteins 33 are responsible for controlling paracellular ion conductance. [34][35][36] As interdigitation increases the cell junction length, it could lead to increased paracellular ion conductance. Indeed, interdigitated cell junctions are observed in leaky epithelia with high paracellular conductance in vivo. Moreover, we provided experimental evidence that the interdigitation may facilitate paracellular ion transport ( Figure S8, supplemental information S6). However, the experiment utilizes blebbistatin, which is not a specific inhibitor of boundary length and may have a broad impact on adherens or tight-junction formation. Further research is required to conclusively demonstrate the regulation of paracellular transport by interdigitation.
). The diffusion coefficient should be proportional to the Young's modulus E and the thickness of the cell boundary L h . d h is also influenced by a compression ratio a.
The mechanical aspect of cell boundaries has been studied in detail. Bambardekar et al. directly manipulate Drosophila cells using an optical tweezer to obtain the mechanical properties like tension and effective viscosity. 37 Cartagena-Rivera et al. measured the MDCK cells using atomic force microscopy (AFM) to measure the tension, effective Young's modulus, and effective viscosity. 38 Diffusion coefficients calculated from these previous studies are 1310 À 6 m 2 =s and 2:8 3 10 À 9 m 2 =s, respectively, which are very different from our estimate (1:0 3 10 À 6 mm 2 =s). This can be understood from the previous analysis that the diffusion coefficient is highly context-dependent-it depends on compression ratio a, which should be small in an interdigitated state.

Stochastic movement of cell boundary
We suggest that the contraction of apical myosin networks focally attached to the junctions may lead to the stochastic movement of cell boundary, but the detail of the movement remains to be elucidated. It is known that pulsed contraction of actomyosin can drive apical constriction (Martin et al., 2009). The dynamic nature of the myosin puncta in the apical cytoplasm implies that its formation may involve myosin-dependent contraction. In agreement with this idea, we observed reduction of the apical myosin network on blebbistatin or Y-27632 treatment. However, it is still not clear whether the puncta push or pull the cell-cell boundary ( Figures 7F-7H). We tried to clarify the relationship between puncta and cell boundaries by mixing wildtype and MRLC2-EGFP cells and observing the cell boundary, but we could observe both pushing and pulling cases. Regardless of whether the puncta push or pull the cell-cell boundary, the loss of the uneven distribution of a-18 by Y-27632 treatment suggests that the uneven tension applied to cell boundaries may be important for the interdigitation formation. In addition, the origin of stochastic movement remains to be elucidated. A model of stochastic pulsatile contraction has been proposed, 39 but they used the complex Swift-Hohenberg equation, 40 which is too abstract and in which the relationship with the molecular mechanism is unclear.

Characteristics of the scaling property of cell-cell boundary
Scaling properties of the growing ragged surface are extensively studied as fractal growth phenomena, and Kardar-Parisi-Zhang (KPZ) equation is frequently utilized to model these phenomena. 41 One difference between these growth phenomena and cell-cell boundary dynamics is that there is no preferential growth direction in the boundary between the same types of cells. Therefore, we used the Edwards-Wilkinson Equation 2, which lacks the nonlinear term of the KPZ equation, which represents the effect of unidirectional growth. Recently, this EW equation was applied to the skull suture 42 and the cell colony boundary. 43 In most cases, the observed scaling exponent is close to À 2, which implies shortening of the boundary length is the plausible cause of scaling. It is known that the scaling exponent in the power spectrum of the steady state of the EW equation is À 2 (supplemental information S2). However, in some cases, we observed that the scaling exponent is considerably different from À 2 ( Figure 5H). One possible explanation is the nonlocal effect of the boundary shape change (supplemental information S3), but further study is needed to elucidate the mechanism of this deviation.

Limitations of the study
Although the ascending limb of the renal tubule has a structure similar to that observed in MDCK cells, we could not conclude that the structure we observed in MDCK cells represents the structure in vivo.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

ACKNOWLEDGMENTS
The authors thank Mami Nakagawa (ABiS) and Yoshimi Yamaguchi (Kyushu University) for technical support, Masayuki Murata (university of Tokyo) for MDCK II cells, Hiroshi Hosoya (Kanagawa University) for pEGFP-N1/MRLC2 plasmid, Akira Nagafuchi (Nara Medical University) for a18 antibodies, and the Developmental Studies Hybridoma Bank for providing rat anti-ZO-1.
We also want to thank Takumi

Lead contact
Further information and requests for resources and reagents should be directed to and will be fullfied by the lead contact, Takashi Miura (miura_t@anat1.med.kyushu-u.ac.jp).

Materials availability
This study did not generate new unique reagents.
Data and code availability d All data reported in this paper will be shared by the lead contact upon request.
d All original code is available in this paper's supplemental information.
d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request. Observation of ex vivo dynamics of renal tubular tight junctions in thick ascending limbs of the mouse Kidneys from young adult R26-ZO-1-EGFP mice (Riken Accession No. CDB0260K 47 ) were cut longitudinally in the plane through the renal pelvis using a scalpel. Half-split kidneys were placed on a glass-bottom dish with the cut side down with DMEM/F12 supplemented with 10% fetal bovine serum. Time-lapse imaging was performed using CellVoyager CV1000 (Yokogawa), a spinning-disk confocal microscope equipped with the stage-top incubator. First, we adjusted the obtained images for brightness and contrast to exclude the effects of fluorescent intensity decrease over time. Next, to remove the effect of tissue-scale deformation, we corrected the cropped time-series images by rigid-body transformation using StackReg plug-in in ImageJ. 55

Numerical simulation
Numerical simulation of the model was implemented by Wolfram Language 12.3 (Wolfram Research Inc.). Partial differential equations of the mathematical models were calculated using explicit Euler scheme. All the source codes were provided as electronic supplemental information.

TER measurement
MDCK II cells were cultured on Transwell filters (0.4 mm pore, polycarbonate; #3401; Corning) for 5 d. 10 mM (-)-Blebbistatin (#021-17041; FUJIFILM Wako) was added to the culture for 4 h, and the electric resistance between the apical and basal chamber was measured using the Millicell ERS-2 electric resistance system (Merck Millipore). The resistance values of the Transwell filters without the cells were measured and subtracted. After subtraction, the unit area resistance was calculated by multiplying the electric resistance by the area of the Transwell filter.

Quantification of the boundary shape
Max projection images of volume data of ZO-1-EGFP MDCK cells were calculated using Fiji. 56 Cell boundary shape was obtained using Ridge Detection plugin. The result was exported as a csv file, and the data was analyzed using Python. 57 The csv file was imported as a dataframe of Pandas, 58 and then, each segment was rotated to standardize as a one-dimensional function hðxÞ using openCV. 59 We manually removed segment data that erroneously contained overhung. Then power spectra of the data were obtained using NumPy, 60 and the results were visulaized using Matplotlib. 61 The source codes were provided as electronic supplemental information. Quantitative data are presented using box-and-whisker plots or histograms using Matplotlib and processed with statistical tests using SciPy unless otherwise stated. Randomization and a priori sample size estimation were not performed. For scaling exponent estimation, linear regression was performed to log-log plots and plotted with 95% confidence interval of scaling exponent, the comparison of exponents was performed using the Mann-Whitney test ( Figures 2C, 5D, 5I, andS4E). For tortuosity and amplitude quantification, Mann-Whitney test was used ( Figures 2D, 2E, 5E, 5F, 5J, 5K, S4F, and S4G). For diffusion coefficient estimation, Kruskal-Wallis test (Figures 8A and 8B) or Mann-Whitney test ( Figure 8C) was used. For MHC-B puncta quantification, (Figure 4C), Mann-Whitney test was used. For TER measurement, data are presented with a bar chart with Student's t-test using Microsoft Excel ( Figure S8B). Details of the statistical analysis are found in the main text.