Mathematical Modeling of Dynamic Cellular Association Patterns in Seminiferous Tubules

Abstract

In vertebrates, sperm is generated in testicular tube-like structures called seminiferous tubules. The differentiation stages of spermatogenesis exhibit a dynamic spatiotemporal wavetrain pattern. There are two types of pattern—the vertical type, which is observed in mice, and the helical type, which is observed in humans. The mechanisms of this pattern difference remain little understood. In the present study, we used a three-species reaction–diffusion model to reproduce the wavetrain pattern observed in vivo. We hypothesized that the wavelength of the pattern in mice was larger than that in humans and undertook numerical simulations. We found complex patterns of helical and vertical pattern frequency, which can be understood by pattern selection using boundary conditions. From these theoretical results, we predicted that a small number of vertical patterns should be present in human seminiferous tubules. We then found vertical patterns in histological sections of human tubules, consistent with the theoretical prediction. Finally, we showed that the previously reported irregularity of the human pattern could be reproduced using two factors: a wider unstable wavenumber range and the irregular geometry of human compared with mouse seminiferous tubules. These results show that mathematical modeling is useful for understanding the pattern dynamics of seminiferous tubules in vivo.

Appendix

1.1 Derivation of the Model

Throughout this manuscript, we used the parameter set which we heuristically chose by coupling the two-component Turing system with a negative feedback loop. However, other parameter sets can reproduce wavetrain patterns. In this appendix, we briefly summarize the derivation of our chosen parameter set and some systematic surveys on possible parameter sets within the framework of the three-component reaction–diffusion Eq. (4). Here, we note that the dispersion relation should satisfy the following conditions:

• 1. $\mathrm {Re}(\lambda _1(0)) < 0$,

• 2. $\mathrm {Re}(\lambda _1(k_{\max })) > 0$,

• 3. $\mathrm {Im}(\lambda _1(k_{\max })) \ne 0$,

• 4. $\mathrm {Re}(\lambda _1(\infty )):=\lim _{k\rightarrow \infty }\mathrm {Re}(\lambda _1(k))<0$.

Condition 1 requires that the spatially uniform mode $k=0$ is stable and does not grow. Conditions 2 and 3 require the existence of the most unstable wavenumber mode with oscillation. Condition 4 requires that the infinite wavenumber mode is stable.

1.1.1 Heuristic Approach

To achieve these conditions, we considered combining the typical two-component Turing system with negative feedback. First, we arbitrarily chose a parameter set $(f_u, f_v,g_u,g_v,D_u,D_v)=(0.5,-1,1.5,-2.1,0.01,1)$ where the two-component reaction–diffusion equation

$$\begin{aligned} \frac{\partial }{\partial t} \left( \begin{array}{c} u \\ v \end{array} \right) = \left( \begin{array}{cc} f_u &{} f_v\\ g_u &{} g_v \end{array}\right) \left( \begin{array}{c} u \\ v \end{array} \right) + \left( \begin{array}{cc} D_u &{} 0\\ 0 &{} D_v \end{array}\right) \frac{\partial ^2}{\partial x^2} \left( \begin{array}{c} u \\ v \end{array} \right) \end{aligned}$$

(12)

shows diffusion-driven instability. Next, we considered another non-diffusive component w which forms a negative feedback loop with component u:

$$\begin{aligned} \frac{\partial }{\partial t} \left( \begin{array}{c} u \\ v \\ w \end{array} \right) = \left( \begin{array}{ccc} f_u &{} f_v &{} \alpha \\ g_u &{} g_v &{} 0\\ -\alpha &{} 0 &{} \beta \end{array}\right) \left( \begin{array}{c} u \\ v \\ w \end{array} \right) + \left( \begin{array}{ccc} D_u &{} 0 &{} 0\\ 0 &{} D_v &{} 0\\ 0 &{} 0 &{} 0 \end{array}\right) \frac{\partial ^2}{\partial x^2} \left( \begin{array}{c} u \\ v \\ w \end{array} \right) . \end{aligned}$$

(13)

Note that the sign of $\alpha$ does not affect the eigenvalues. In the limit of $k\rightarrow \infty$, the characteristic polynomial of the system (13) should be

$$\begin{aligned} k^4 D_u D_v (\beta -\lambda ) + (\mathrm {l.o.t.\ of\ }k) = 0. \end{aligned}$$

Thus, to satisfy condition 4, $\mathrm {Re}(\lambda _1( \infty )) = \beta < 0$. We numerically searched $0\le \alpha \le 10$ and $-10\le \beta \le 0$, and found the possible parameter range (Fig. 12a). We chose $(\alpha ,\beta )=(1,-0.1)$ in the possible parameter range. We also considered that w forms a negative feedback loop with component v, not u, but could not find $(\alpha ,\beta )$ satisfying conditions 1–4 in the same range.

Fig. 12

The derivation and systematic surveys of model parameters in three-component reaction–diffusion equations for wavetrain pattern formation. a The possible parameter range for the negative feedback parameter $\alpha$ and the autoregulation of a non-diffusion component $\beta$ in the heuristic approach. The possible range is shown in black. b, c Random search of reaction term parameters between -1 and 1 when $(D_u, D_v, D_w) = (0.01, 1, 0)$. b The histogram for each parameter. c The sign similarity index for each parameter pair. Red and blue correspond to the same and opposite signs, respectively. d, e Full search of reaction term parameters when $(D_u, D_v, D_w) = (0.01, 0.01, 1)$. d The histogram for each parameter. e The sign similarity index for each parameter pair

1.1.2 Random Search of Reaction Term Components with the Fixed Diffusion Coefficient Set

Based on the parameter sets found by the heuristic approach, we surveyed the possible parameter range for all the coefficients in the linear part of the reaction term in a three-component reaction–diffusion equation. We employed the same diffusion coefficients as the heuristic approach, $(D_u, D_v, D_w) = (0.01,1,0)$. Each component of the reaction term matrix A was randomly chosen from the uniform distribution $(-1,1)$.

We tried 2,000,000 random parameter sets and 1,201 sets were found to satisfy all conditions (0.06%). In the nine parameters, $f_u$ and $g_v$ were biased to positive and negative values, respectively (Fig. 12b). It seems u and v form the Turing system since $D_u < D_v$. However, only 27% (327 sets) satisfied the conditions required for Turing instability ($f_u + g_v < 0$, $f_u g_v - f_v g_u > 0$, $D_v f_u + D_u g_v > 0$, $(D_v f_u + D_u g_v)^2 - 4 D_u D_v (f_u g_v - f_v g_u) > 0$). Also, $h_w$ was always negative, the same as in the heuristic approach (Fig. 12b).

We also assessed the combination of two parameters, by introducing the sign similarity index as the frequency of a particular parameter pair taking the same sign (Fig. 12c). Notably, $f_w$ and $h_u$ always showed the opposite signs. These connections would form a negative feedback loop between u and w, consistent with our heuristic approach.

1.1.3 Full Search

We also surveyed other diffusion coefficient sets. Because of the high-dimensional parametric space, we limited diffusion coefficients to either 0, 0.01, or 1, and reaction matrix components to either $-1, -0.1$, 0, 0.1, or 1. We excluded cases where $D_u=D_v=D_w$ and only considered $D_u \le D_v \le D_w$. Only 876 of 13,671,875 cases showed the desired dispersion relation (0.006%). Most of them were either $(D_u, D_v, D_w) = (0, 0.01, 1)$ (284 cases) or $(D_u, D_v, D_w) = (0.01, 0.01, 1)$ (424 cases). The analyses in the previous section included the former cases. For $(D_u, D_v, D_w) = (0.01, 0.01, 1)$, we also surveyed the reaction term parameters (Fig. 12d, e). Only diagonal components $f_u$, $g_v$, and $h_w$ were biased to either positive or negative. Similar to the previous section, the autoregulation of the fastest diffusion component (w in this case) tended to be negative (Fig. 12d). In contrast, those of the slow diffusion components u and v were slightly biased to be positive (Fig. 12d). We also identified an important negative feedback loop between two slow diffusion components ($f_v$ and $g_u$, Fig. 12e). However, neither u and w, or v and w constituted a Turing system for any cases.

1.2 Numerical Simulations with an Alternative Parameter Set

To clarify that our wavetrain pattern arguments are not largely determined by our specific parameter set, we performed a set of numerical simulations using the three-species reaction–diffusion model (3) where

$$\begin{aligned} A = \left( \begin{array}{ccc} 2.9 &{} -2.9 &{} -2 \\ 1.5 &{} -1 &{} 0 \\ 100 &{} 0 &{} -100 \end{array}\right) , D = \left( \begin{array}{ccc} 0.024 &{} 0 &{} 0 \\ 0 &{} 0.006 &{} 0\\ 0 &{} 0 &{} 3 \end{array} \right) , C(\mathbf {u}) = \left( \begin{array}{ccc} u^3 \\ 0\\ 0 \end{array} \right) . \end{aligned}$$

(14)

This parameter set satisfies conditions shown in Fig. 2.

We considered the two-dimensional model with the same domain size as in the main text ($25.6 \times 3.2$). In this simulation, the pattern size is shorter than the circumferential length as $d = 1.71 < L_0$. Thus, we expected both vertical and helical patterns would appear ((i) and (iii) in Fig. 7). We confirmed both patterns appeared with the numerical simulations (Fig. 13).

Fig. 13

Numerical simulations of a two-dimensional model (9) with an alternative parameter set (14). Consistent with the theoretical predictions, vertical and helical patterns appeared. Mixed and segmented patterns were also observed

1.3 Weakly Nonlinear Analysis of a One-Dimensional Model

We examined the dynamics of the amplitudes of forward and backward waves by weakly nonlinear analysis. First, we defined $\mathbf {A_b}$ and $\mathbf {A_f}$, which represent the backward and forward waves, respectively:

$$\begin{aligned} \begin{aligned} \mathbf {A_b}&= \varvec{\omega }_1 \mathrm{e}^{i k(x + c t)} + \varvec{\omega }_2 \mathrm{e}^{-i k (x + c t)}\\ \mathbf {A_f}&= \varvec{\omega }_1 \mathrm{e}^{-i k(x - c t)} + \varvec{\omega }_2 \mathrm{e}^{i k (x - c t)}. \end{aligned} \end{aligned}$$

(15)

We also defined $a_\mathrm{b}$ and $a_\mathrm{f}$ as amplitudes of the backward and forward waves, respectively.

$$\begin{aligned} \begin{aligned} \mathbf {u}&= a_\mathrm{b} \mathbf {A_b} + a_\mathrm{f} \mathbf {A_f}\\&=\varvec{\omega }_1(a_\mathrm{b} \mathrm{e}^{ik(x+ct)} + a_\mathrm{f} \mathrm{e}^{-ik(x-ct)}) + \varvec{\omega }_2(a_\mathrm{f} \mathrm{e}^{ik(x-ct)} + a_\mathrm{b} \mathrm{e}^{-ik(x+ct)})\\&= (\varvec{\omega }_1 \mathrm{e}^{i k(x+ct)} + \varvec{\omega }_2 \mathrm{e}^{-i k (x+ct)})a_\mathrm{b} + (\varvec{\omega }_1 \mathrm{e}^{-ik(x-ct)} + \varvec{\omega }_2 \mathrm{e}^{ik(x-ct)})a_\mathrm{f}. \end{aligned} \end{aligned}$$

(16)

The left-hand side of the model (3) was:

$$\begin{aligned} \begin{aligned}&ikc(\varvec{\omega }_1 \mathrm{e}^{i k(x+ct)} - \varvec{\omega }_2 \mathrm{e}^{-i k (x+ct)})a_\mathrm{b}+(\varvec{\omega }_1 \mathrm{e}^{i k(x+ct)} + \varvec{\omega }_2 \mathrm{e}^{-i k (x+ct)})\frac{\partial a_\mathrm{b}}{\partial t}\\&\quad +ikc(\varvec{\omega }_1 \mathrm{e}^{-i k(x-ct)} - \varvec{\omega }_2 \mathrm{e}^{i k (x-ct)})a_\mathrm{f}+(\varvec{\omega }_1 \mathrm{e}^{-ik(x-ct)} + \varvec{\omega }_2 \mathrm{e}^{ik(x-ct)})\frac{\partial a_\mathrm{f}}{\partial t}. \end{aligned} \end{aligned}$$

(17)

The right-hand side of the model (3) was:

$$\begin{aligned} \begin{aligned}&(A-k^2 D) \mathbf {u} - C(\mathbf {u})\\&\quad = \lambda _1 \varvec{\omega }_1(a_\mathrm{b} \mathrm{e}^{ik(x+ct)} + a_\mathrm{f} \mathrm{e}^{-ik(x-ct)})\\&\qquad +\lambda _2 \varvec{\omega }_2(a_\mathrm{f} \mathrm{e}^{ik(x-ct)} + a_\mathrm{b} \mathrm{e}^{-ik(x+ct)})-C(\mathbf {u}). \end{aligned} \end{aligned}$$

(18)

The nonlinear term $C(\mathbf {u})$ was expanded as follows

$$\begin{aligned} \begin{aligned} C(\mathbf {u})&= C(a_\mathrm{b} \mathbf {A_b} + a_\mathrm{f} \mathbf {A_f})\\&=a_\mathrm{b}^3 \mathbf {A_b}^3 + 3 a_\mathrm{b}^2 a_\mathrm{f} \mathbf {A_b}^2 \mathbf {A_f} + 3 a_\mathrm{b} a_\mathrm{f}^2 \mathbf {A_b} \mathbf {A_f}^2 + a_\mathrm{f}^3 \mathbf {A_f}^3. \end{aligned} \end{aligned}$$

(19)

Using (18) and (19), we obtained the right-hand side of the model (3):

$$\begin{aligned}&\lambda _1 \varvec{\omega }_1(a_\mathrm{b} \mathrm{e}^{ik(x+ct)} + a_\mathrm{f} \mathrm{e}^{-ik(x-ct)}) +\lambda _2 \varvec{\omega }_2(a_\mathrm{f} \mathrm{e}^{ik(x-ct)} + a_\mathrm{b} \mathrm{e}^{-ik(x + ct)})\\&-a_\mathrm{b}^3 \mathbf {A_b}^3 - 3 a_\mathrm{b}^2 a_\mathrm{f} \mathbf {A_b}^2 \mathbf {A_f} - 3 a_\mathrm{b} a_\mathrm{f}^2 \mathbf {A_b} \mathbf {A_f}^2-a_\mathrm{f}^3 \mathbf {A_f}^3. \end{aligned}$$

Substituting $\mathbf {A_f}$ and $\mathbf {A_b}$ for (6), we obtained:

$$\begin{aligned} \begin{aligned}&\lambda _1 \varvec{\omega }_1(a_\mathrm{b} \mathrm{e}^{ik(x + ct)} + a_\mathrm{f} \mathrm{e}^{-ik(x-ct)}) +\lambda _2 \varvec{\omega }_2(a_\mathrm{f} \mathrm{e}^{ik(x-ct)} + a_\mathrm{b} \mathrm{e}^{-ik(x+ct)})\\&\quad -((3 a_\mathrm{b}^3 \varvec{\omega }_1^2 \varvec{\omega }_2 + 6 a_\mathrm{b} a_\mathrm{f}^2 \varvec{\omega }_1^2 \varvec{\omega }_2) \mathrm{e}^{ik(x+ct)} +(3 a_\mathrm{f}^3 \varvec{\omega }_1^2 \varvec{\omega }_2 + 6 a_\mathrm{b}^2 a_\mathrm{f} \varvec{\omega }_1^2 \varvec{\omega }_2) \mathrm{e}^{-ik(x-ct)}\\&\quad +(3 a_\mathrm{f}^3 \varvec{\omega }_1 \varvec{\omega }_2^2 + 6 a_\mathrm{b}^2 a_\mathrm{f} \varvec{\omega }_1 \varvec{\omega }_2^2) \mathrm{e}^{ ik(x-ct)} +(3 a_\mathrm{b}^3 \varvec{\omega }_1 \varvec{\omega }_2^2 + 6 a_\mathrm{b} a_\mathrm{f}^2 \varvec{\omega }_1 \varvec{\omega }_2^2) \mathrm{e}^{-ik(x + ct)})\\&\quad +(\text {h.o.t.}). \end{aligned} \end{aligned}$$

(20)

We defined a real vector as: $\varvec{\Omega } = \varvec{\omega }_1 \varvec{\omega }_2$, which represented the squared absolute value of components $\varvec{\omega }_1$ and $\varvec{\omega }_2$. Then, (20) was simplified as follows:

$$\begin{aligned} \begin{aligned}&\lambda _1 \varvec{\omega }_1(a_\mathrm{b} \mathrm{e}^{ik(x + ct)} + a_\mathrm{f} \mathrm{e}^{-ik(x-ct)}) +\lambda _2 \varvec{\omega }_2(a_\mathrm{f} \mathrm{e}^{ik(x-ct)} + a_\mathrm{b} \mathrm{e}^{-ik(x+ct)})\\&\quad -( \varvec{\Omega }\varvec{\omega }_1 ((3 a_\mathrm{b}^3 + 6 a_\mathrm{b} a_\mathrm{f}^2) \mathrm{e}^{ik(x+ct)}+ (3 a_\mathrm{f}^3 + 6 a_\mathrm{b}^2 a_\mathrm{f}) \mathrm{e}^{-ik(x-ct)})\\&\quad + \varvec{\Omega }\varvec{\omega }_2( (3 a_\mathrm{f}^3 + 6 a_\mathrm{b}^2 a_\mathrm{f}) \mathrm{e}^{ik(x-ct)}+ (3 a_\mathrm{b}^3 + 6 a_\mathrm{b} a_\mathrm{f}^2) \mathrm{e}^{-ik(x + ct)}))\\&\quad +(\text {h.o.t.}). \end{aligned} \end{aligned}$$

(21)

To assay the contribution of a nonlinear term ($\varvec{\Omega }\varvec{\omega }_1, \varvec{\Omega }\varvec{\omega }_2$) to each eigenvector component, we defined $\alpha _1$, $\alpha _2$, $\alpha _3$ as follows:

$$\begin{aligned} \varvec{\Omega }\varvec{\omega }_1 = \alpha _1 \varvec{\omega }_1 + \alpha _2 \varvec{\omega }_2 + \alpha _3 \varvec{\omega }_3. \end{aligned}$$

(22)

Using the characteristics of complex conjugates, $\varvec{\Omega }\varvec{\omega }_2$ was defined as follows:

$$\begin{aligned} \varvec{\Omega }\varvec{\omega }_2 = \varvec{\Omega }\bar{\varvec{\omega }}_1 = {\bar{\alpha }}_1 \bar{ \varvec{\omega }}_1 + {\bar{\alpha }}_2 \bar{\varvec{\omega }}_2 + {\bar{\alpha }}_3 \bar{\varvec{\omega }}_3. \end{aligned}$$

(23)

We considered only $\varvec{\omega }_1 \mathrm{e}^{-ik(x-ct)}$, $\varvec{\omega }_2 \mathrm{e}^{ik(x-ct)}$, $\varvec{\omega }_1 \mathrm{e}^{-ik(x+ct)}$ and $\varvec{\omega }_2 \mathrm{e}^{-ik(x+ct)}$, as all other components will be diminished because of the effect of the linear term. In addition, we focused on the forward wave components $\mathrm{e}^{-ik(x-ct)}$ and $\mathrm{e}^{ik(x-ct)}$ because of the symmetry.

By calculating (17) $=$ (18) using $kc=\lambda _i$, we obtained:

$$\begin{aligned}&\frac{\mathrm{d}}{\mathrm{d}t} a_\mathrm{f} = a_\mathrm{f}\left( \lambda _r - \alpha _r\left( 3 a_\mathrm{f}^2 + 6 a_\mathrm{b}^2\right) \right) \end{aligned}$$

(24)

where $\lambda _r$ and $\alpha _r$ are the real parts of $\lambda _1$ and $\alpha _1$, respectively. It is, actually, a natural consequence from the system’s reverse symmetry that the ratio between the internal and external interactions in (24) is 3 to 6.

1.4 Effect of Finite Domain Size on Pattern Selection

We determined the effect of finite domain size upon pattern selection. Because domain size is limited in numerical simulations, there are discrete numbers of possible solutions, which make pattern selection more complex. The patterns were classified into at least five categories, as shown in Fig. 6c. These patterns reflected the effect of finite domain size along the long axis.

${\mathbf{P}}_1$: When $s \ge 1.4$, there were only vertical patterns (i). The pattern was restricted by the relationship between the pattern size and circumferential length in this range.

$\mathbf{P}_2$: When $s=1.2$, the horizontal pattern (ii) appeared. In this range, horizontal patterns were formed by the relationship between the pattern size and circumferential length, and were the most unstable compared with vertical or helical patterns (orange point in Fig. 14a).

$\mathbf{P}_3$: When $0.7 \le s \le 1$, helical patterns were highly dominant. For instance, when $s = 0.9$, there are two wavelengths allowed in helical waves (gray points in Fig. 14b and one wavelength allowed in horizontal waves (blue point in Fig. 14b) near the most unstable wavelength. Therefore, we expected that the helical pattern was much more frequent in this parameter range.

$\mathbf{P}_4$: When $0.3 \le s \le 0.6$, the number of helical patterns was as almost twice that of vertical patterns. For instance, when $s = 0.4$, one gray point and one blue point represented the wavelengths of helical and vertical waves near the most unstable wavelength, respectively (Fig. 14c, arrows). Therefore, the left helical, right helical, and vertical patterns appeared at the approximately the same frequency in this parameter range.

$\mathbf{P}_5$: When $s \le 0.2$, helical patterns were highly dominant. In this range, helical patterns in which one cycle of differentiation corresponded to the circumferential length (Fig. 7a(iii)), and also helical patterns in which two or more modes corresponded to the circumferential length (Fig. 7a(iv)) were allowed.

Fig. 14

(Color figure online) Dispersion relation and possible solutions when $s=1.2$ (a), $s = 0.9$ (b), and $s = 0.4$ (c) with the finite domain size. Blue, orange, and gray points represent the possible wavelengths of vertical, horizontal, and helical waves, respectively. Arrows mark the points near $k_{\max }$