Modeling the SARS-CoV-2 Omicron variant dynamics in the United States with booster dose vaccination and waning immunity

1.   Introduction

Pancreatic development presents fascinating topic of how cell fate decisions are realized. It inspires many scholars to explore relevant topics from both experimental and theoretical aspects. There are three cell types: exocrine cells, endocrine cells, and ductal cells in the mature pancreas organ which derived from the embryonic pancreas ^[1,2,3]. The vertebrate pancreas comprises of two main components: the endocrine and exocrine compartment. The endocrine cell secretes hormones in the bloodstream and controls glucose homeostasis. However, the exocrine compartment secretes the main digestive enzymes in the gut lumen ^[4]. In addition, the islets of Langerhans in endocrine tissue, further develop into four pancreatic endocrine cell types: $\alpha$-, $\beta$-, $\delta$-, and $PP$-cells which secret glucagon, insulin, somatostatin, and pancreatic polypeptide, respectively ^[4,5]. The exocrine tissue is composed of acini and ducts which specialized in enzymes production and secreting bicarbonated water, vehiculating the pancreatic enzymes to the intestine, respectively ^[4]. As the lack of insulin may cause diabetes and pancreatic cancer, much attention has been focused on the regulatory mechanisms about the differentiation of each cell type ^[6,7].

One of the early process of pancreatic cell development is a fate decision between endocrine and exocrine, in which a number of proteins proteins are involved during the development of the mouse pancreas, the protein Hnf6 is expressed in the epithelial cells. And these epithelial cells represent precursors cells of the exocrine and endocrine pancreatic. The pancreatic endocrine development at the precursor process is regulated by Hnf6 and it is identified as the first positive regulator of the pro-endocrine protein neurogenin3 (Ngn3) in the pancreas ^[8]. While the Ngn3 is expressed in cells which are destined to become endocrine cells, the islets of Langerhans ^[9]. Ngn3 has an ectopic expression lead to cells to an endocrine fate in the early pancreatic bud, but the almost completely $\alpha$-cells are obtained^[10,11] and additional signals are need to deflect these cells into alternate fates such as $\beta$-cells. Previous findings suggest that Isl1 expression is need to the development of islet cells in pancreatic epithelial cells and it is a terminal endocrine fate marker downstream of Ngn3^[9,12]. $Ptf1a$ is one of the gene necessary for exocrine differentiation and it was was first identified in the nucleus of exocrine cultured pancreatic cells in an extensive footprint analysis of the of digestive enzymes promoters^[13] and its isolation can obtain the genetic characteristics encoding basic helix-loop-helix protein of 48 kDa ^[14]. Ptf1a is commonly used as a marker of both the precursor cells and differentiated status exocrine cells of pancreatic in zebrafish ^[4]. It is the fate specific marker of exocrine progenitor cells ^[15]. In addition, PTF1L-complex is necessary for acinar cell differentiation, i.e., RBPJ and Ptf1a bind to a small basic helix-loop-helix protein to form a trimeric PTF1J-complex which also activates the expression of Rbpjl and leads to an increase in protein that progressively instead of RBPJ in the trimeric PTF1J-complex, forming the PTF1L-complex ^[16].

Notch signaling is necessary for many key regulatory events in cellular destiny during cell development, such as nervous system, skin and so on ^[17]. In addition, previous studies indicate that Notch signaling also plays a key roles in the development of the pancreas and is important in determining the endocrine and progenitor/exocrine destiny of developing pancreas ^[18,19]. The reduction of Notch can trigger an increase in expression of the pro-endocrine gene $Ngn3$ and lead to differentiation toward endocrine destiny direction. However, at the normal level of Notch, cells express Hes1 and p48 and select the direction of exocrine destiny ^[18]. Related study also indicates that Notch signalling activates Ptf1 through the direct interaction between RBP-Jk, NotchIC and p48, which leads to the selection towards the exocrine destiny direction ^[19]. Moreover, endocrine progenitors are generated by lateral specification through Notch signalling involving Delta/Notch (D/N) during pancreas development^[20]. The progenitor-like cells (adult acinar cells) convert into $\beta$-cells upon cell-cell contacts disruption through enzymatic tissue separation which involves inactivation of Notch signaling ^[21,22].

Notch receptor exists on the cell membrane of a certain cell binds to a ligand Delta-like1 (Dll1) on the surface of other contact cells to produce the Notch intracellular domain (NICD) ^[23]. It is known that Notch exists on the cell membrane can be induced by ligand (Dll1) which comes from in neighboring cells (trans-activation) and also be suppressed by Dll1 in the same cell (cis-inhibition) ^[24]. The cis-inhibition usually affects the direction of D/N signaling and should be incorporated into specification models of lateral suppression^[25]. In addition, trans-activation has received extensive attention both experimentally and theoretically, for example, related studies in the control of vertebrate neurogenesis ^[26] and somite formation ^[27]. The pancreatic cell fate decision ^[18] and vertebrate neurogenesis ^[26] usually involve Notch signaling mechanism which is called lateral suppression. A unified lateral suppression mode of D/N regulated neurogenic and pancreatic progenitor specification has been mentioned ^[20]. Therefore, it is an important and meaningful topic to analyze fate decisions of pancreatic cells based on Notch signaling.

Mathematical modeling and analysis can help us understand regulatory and developmental mechanism in cell fate decisions ^[28,29]. Some studies explored the mechanism of pancreatic cell fate decision via mathematical models in recent years, e.g., predicting pancreas reprogramming with a stratification multiple attractor model ^[15], understanding transdifferentiation of pancreatic cells with a minimal model and contact-mediated signaling pathways ^[22], a dynamical model of a growing duct was established, which leads to an oscillatory phase before the decision of endocrine progenitors fate by lateral suppression ^[20]. Here, we construct mathematical models at the single cell level and multicellular level for pancreas destiny choice which involves dynamics of winding with trans-activation and cis-inhibition within Notch signaling. We further analyze the models by means of bifurcation analysis and understand the regulatory mechanism of pancreas fate decisions.

2.   Methods

Gene regulatory networks can be mathematically modeled and analyzed in terms of differential equations. It can help to understand the complex feedback mechanisms underlying cell fate decision. We construct a minimal regulatory network, as shown in Figure 1(A), to study pancreatic cell fate decisions by combining related literatures, as shown in Table 1. The expression of Hnf6 is suppressed by Ptf1a in mature exocrine cells (Ptf1a expression is a hallmark of mature exocrine cells) ^[15], and as a mature exocrine cell, i.e., mature acinar cells (Ptf1a and Rbpjl expression are hallmarks of mature acinar cells) ^[30]. So we have a possible inference that Rbpjl inhibits Hnf6, and the reference ^[56] further confirm this inference.

The state of each cell type is specified by five key proteins: Hnf6, Ngn3, Ptf1a, Isl1, and Rbpjl. More specifically, epithelial cells are commonly regarded as precursors of the exocrine and endocrine cells and Hnf6 is expressed in these cells. It is known that the upstream protein Hnf6 either directly or indirectly activates Ngn3 and Ptf1a ^[22], and under normal circumstances, the adult pancreas does not express these proteins (Ngn3 and Ptf1a) at an advanced stage of development ^[9]. Ngn3 is a marker of pro-endocrine protein which is transiently existence in early pancreatic development. When Ngn3 is induced in the progenitor cell within the pancreas, the cell is destined to be an endocrine (islet) cell ^[5]. In addition, the expression of Isl1 in pancreatic epithelial cells is essential for the development of islet cells ^[9,12]. Isl1 is generally regarded as a terminal endocrine destiny marker downstream of Ngn3, and once activated, keeps its expression by positive self-activation ^[22]. Isl1 as an islet cell maturation protein, its role is to inhibit the expression of upstream protein Hnf6. Ptf1a is the only known protein that is indispensable and sufficient to active the exocrine cell destiny ^[4,14], and it is the fate specific marker of exocrine progenitor cells ^[15]. In mature acinar cells, Ptf1a establishes an auto-regulation via the trimeric complex PTF1 to reinforce and maintain its own expression. Indeed, PTF1 forms dual self-regulating loops of Ptf1a and Rbpjl which may maintain the stability of pancreatic acinar cell phenotypes ^[30]. Moreover, Rbpjl represents a terminal exocrine fate marker, and it is necessary for acinar cells to be fully mature^[31].

Furthermore, to describe and understand the development of progenitor cells into endocrine and exocrine cells easily, we divide the differentiation process into two stages, as shown in Figure 1(B). At the first stage, the lineage specification of pancreatic progenitor cells into endocrine or exocrine progenitor cells occurs. While at the second stage, the lineage specification of endocrine and exocrine progenitor cells into mature islet and acinar cells occurs, respectively. The mathematical model for the network shown in Figure 1(A) is described by Eq (1). The definitions of variables and parameters in Eq (1) are given in Table 3 in the Appendix.

Moreover, to reveal the mechanism of Notch signaling in pancreatic cell fate decision, we construct a multicellular mathematical model by combining the minimal model Eq (1) and Notch signaling, as shown Figure 2. Mathematical models involving Notch signaling have been presented in ^[24,48]. The model here involves the following aspects: First, Notch in cell-$i$, i.e., $x_{7, i}$, combined with extracellular Dll1 with concentration $\langle x_{6, j}\rangle_i$, i.e., the average Dll1 level over all adjacent cells $j$ of $i$, resulting in the generation of the NCID, $x_{8, i}$. In the same manner, Notch in adjacent cells, $\langle x_{7, j}\rangle_i$ can combined with Dll1. And, Notch can also combined with Dll1 within the same cell, resulting in inactivation of Notch. The multicellular system can be expressed by Eq (2) for the concentrations of free Notch, $x_{7, i}$, free Dll1, $x_{6, i}$, NICD, $x_{8, i}$, Hnf6, Ngn3, Ptf1a, Isl1, and Rbpjl in cell $i$ $(i = 1, ..., N, \ j\neq i)$. The definitions and values of all parameters in Eq (2) are given in Table 4 in the Appendix.

3.   Results

3.1.   Regulatory mechanism of pancreatic cell fate decision

In this section, we mainly reveal the mechanism of pancreatic cell fate decision through Eq (1), e.g., pancreatic progenitor cell differentiation, endocrine or exocrine progenitor cell differentiation and transdifferentiation of acinar cell. In addition, We also explore the role of auto-regulation loops in pancreatic cell fate decision, and the results show that the autoregulatory loops play a role in maintaining the stability of mature acinar cells.

3.1.1.   Decision switches during pancreatic cell differentiation

The process of pancreatic cell differentiation can be explored and understood by the dynamics of Eq (1). In fact, pancreatic progenitor cells first develop into transitional cells, e.g., endocrine and exocrine progenitor cells, and finally differentiate into mature islet and acinar cells. Now, the mechanism of decision switches is explored by bifurcation analysis of Eq (1). At the first stage, the pancreatic progenitor lineage branches into exocrine and endocrine progenitor lineages, and the binary decision switch is induced by the mutually suppressing protein pair Ptf1a and Ngn3 ^[15,33]. The bifurcation diagram with $a_1 = a_2 = a_{NP}$ as the control parameter is shown in Figure 3(A). At the critical value $a_{NP} = a_{c}$, a pitchfork bifurcation occurs. When $a_{NP} < a_{c}$, cells are at the pancreatic progenitor lineage. As $a_{NP}$ increases, i.e., when $a_{NP} > a_{c}$, the progenitor cells develop into exocrine or endocrine progenitor cells. Cells with high Ngn3 and low Ptf1a expression are endocrine progenitor cells, i.e., $S_{A}$. While cells with low Ngn3 and high Ptf1a are exocrine progenitor cells, i.e., $S_{B}$, as shown in Figure 3(B). The decision switch is realized through occurrence of the pitchfork bifurcation.

At the second stage, endocrine or exocrine progenitor cells can further transform into adult islet or acinar cells, respectively. The fate decision switch is realized through saddle-node bifurcation which is produced by the coupling between positive auto-regulation of Isl1 and negative loop between Ngn3 and Isl1, or the coupling between positive auto-regulation of Rbpjl and negative loop between Ptf1a and Rbpjl, respectively. Pancreatic progenitor cells transform into either endocrine or exocrine progenitor cells after the first stage. If pancreatic progenitor cells choose to differentiate into endocrine progenitor cells, the ultimate fate will be islet cells. Mathematically, when $a_4 < a_{4_c}$, the system stays at the endocrine progenitor state. As $a_4$ increases, i.e., when $a_4 > a_{4_c}$, endocrine progenitor cells transform into adult islet cells, as shown in Figure 4. If pancreatic progenitor cells choose to differentiate into exocrine progenitor cells, then the ultimate fate will be acinar cells, i.e., when $a_6 < a_{6_c}$, they stays at the exocrine progenitor state, as $a_6$ increases, i.e., when $a_6 > a_{6_c}$, exocrine progenitor cells turn into adult acinar cells, as shown in Figure 5.

Remark 1. We simulate the transformation of pancreatic progenitor cells into endocrine or exocrine progenitor cells by small perturbing initial values of Eq (1). Such a kind of perturbation is performed to mimic transient induction of differentiation signals. When they choose to differentiate into endocrine progenitor cells with high Ngn3 and low Ptf1a, i.e., $S_A$ in Figure 3(B), the ultimate fate will be islet cells. When they choose to differentiate into exocrine progenitor cells with high Ptf1a and low Ngn3, i.e., $S_B$ in Figure 3(B), due to the autoregulation of Ptf1a, the high level of Ptf1a ($S_B$) induces much higher Ptf1a (its level near to 2 in Figures 5(A) and 6(A)), then the ultimate fate will be acinar cells. Between the first and second stages, the state with high Ptf1a or much higher Ptf1a corresponds to exocrine progenitor cells, as mentioned in Table 2.

3.1.2.   Auto-regulation of Ptf1a and Rbpjl is necessary

Auto-regulation loops are common in regulatory networks which may control the behavior and maintenance of stability in biological systems ^[34,35]. In development cells, such a simple regulatory pattern acts as molecular memory to keep the mature phenotype for specific cell lineage, e.g., auto-regulations of Ptf1a and Rbpjl keep the phenotype of mature acinar cells ^{[30,31,32,36,37]}. Now, we reveal the mechanism of acinar cell fate decision by bifurcation analysis of Eq (1). For convenience, we set the auto-activation strength of Ptf1a and Rbpjl to be the same, i.e., $a_3 = a_7 = a_P$. Bifurcation diagrams of Eq (1) with $a_P$ as the control parameter is shown in Figure 6. When $a_P$ is small, i.e., $a_P < a_{P_{c_1}}$, the system stays at the exocrine progenitor cell state. When $a_{P_{c_2}} < a_P < a_{P_{c_3}}$, Ptf1a is superinduced with higher expression, but the expression of Rbpjl is still low. The superinduction of Ptf1a expression is necessary and important developmental event for the formation of mature acinar cells ^[30]. As $a_P$ increases further, i.e., when $a_P > a_{P_{c_3}}$, the exocrine progenitor cells eventually transform into mature acinar cells and maintain the high expression of Ptf1a and Rbpjl, as shown in Figure 6.

In fact, the protein Ptf1a binds to a common E-proteins (e.g., TCF12, HEB) and RBPJ to produce a trimeric complex PTF1-J ^[16,30,36], which binds to the $5^{'}$ enhancer to auto-activate Ptf1a. The enhancer and Ptf1a form an auto-regulatory loop to enhance and keep the expression of Ptf1a ^[30]. Along with the $5^{'}$ enhancer, Ptf1a obtains superinduction via the activation of unknown activators, and it may involve increased auto-activation strength of Ptf1a, i.e., $a_3$, leading to its superinduced expression ($a_{P_{c_2}} < a_P < a_{P_{c_3}}$ in Figure 6(A)). When the Ptf1a is elevated, the protein Rbpjl is activated by the PTF1-J complex. The RBPJL (protein translated by Rbpjl) substitutes RBPJ within the PTF1 complex to produce PTF1-L, which binds to Ptf1a enhancer and the Rbpjl promoter, and keeps their activity, inducing an increase in auto-activation of Rbpjl, i.e., $a_7$. Indeed, the trimeric complex PTF1-L constitutes dual auto-regulatory loops to keep the expression of Ptf1a and Rbpjl ($a_P > a_{P_{c_3}}$ in Figure 6) ^[30]. Since the auto-regulations of Ptf1a and Rbpjl are regulated by the PTF1 complex (PTF1-J, PTF1-L), for convenience, we set the auto-activation strength of Ptf1a and Rbpjl to be the same, i.e., $a_P$.

3.1.3.   Mechanism of acinar-to-islet cell conversion

During embryonic evolution, cells gradually turn into more professionalization. But, many studies clearly show that adult differentiated cells still have the ability to change the fate of cells under appropriate conditions ^[38,39]. New technology in regenerative medicine aim to take advantage of the cell plasticity to take the place of diseased tissue by directional transformation of cells from cells of other types. Transdifferentiation, also known as the indirect change in cell specification from one cell lineage to another, often involves a phase of dedifferentiation to restore multipotency. However, it is also possible to drive cells to change lineages directly ^[40]. Acinar cells are an candidate source of transdifferentiation due to the common developmental source of exocrine and endocrine cells and the large number of acinar cells in pancreas. To investigate the mechanism of acinar-to-islet cell transformation, we study the stability of the system and their dependence on parameter values via executing bifurcation analysis. Figure 7 shows that Isl1 and Rbpjl have three stable states within a vast range of parameter values. The three cell types are acinar cell type (high Rbpjl, low Isl1), islet cell type (low Rbpjl, high Isl1), and progenitor-like multipotent type (the intermediate solid blue line). When the initial state is at the mature acinar type at $a_P > a_{P_{s_1}}$, as $a_P$ decreases, the acinar dedifferentiation towards a progenitor-like multipotent (exocrine progenitor) cell type occurs within $a_{P_{s_2}} < a_P < a_{P_{s_1}}$. At $a_P < a_{P_{s_2}}$, the progenitor-like multipotent cells redifferentiate into islet cells. In fact, the trimeric complex (PTF1-L) forms dual auto-regulatory loops, and when the complex is inhibited via genetic manipulation, auto-activation is weakened due to the decrease of the auto-activation strength $a_P$. In other words, a small $a_P$ value indicates a large degree of damage to the auto-regulatory loops of ptf1a and Rbpjl. Therefore, below the critical value $a_{P_{s_2}}$, only the islet cell fate keeps stable.

Acinar cells are an candidate source of transdifferentiation due to the co-developmental source of exocrine and endocrine cells and the large number of acinar cells in pancreas ^[22]. Indeed, the transdifferentiation of acinar cells into new $\beta$-cells was confirmed in vivo experiment in mice by manipulating ectopic expression of pivotal transcription factors ^[41]. New $\beta$-cells could be obtained via transdifferentiation of acinar cells into islet cells because islet cells could differentiates into one of endocrine cell subtypes, i.e., $\beta$-cells. Excitingly, such transdifferentiation has also been confirmed in vitro cultures without genetic operation, i.e., using only microenvironmental changes, e.g., enzymatic tissue dissociation ^[42,43]. The transdifferentiation via key transcription factor manipulation and the microenvironmental changes may be involved in preventing the formation of the auto-regulatory loops of ptf1a and Rbpjl so that acinar cells cannot maintain the expression of Ptf1a, thus undergoing transdifferentiation to islet cells and eventually differentiating into $\beta$-cells. The explanation of how enzymatic tissue dissociation prevents the formation of the auto-regulatory loops of ptf1a and Rbpjl can be as follows. The downstream protein connexin of Mist (mentioned in Section {3.1.4}) is important for the formation of the auto-regulatory loop of Ptf1a and when the connexin structure is broken via enzymatic tissue dissociation, the auto-regulatory loops of ptf1a and Rbpjl are subsequently impaired.

3.1.4.   The PTF1 complex and Mist1 may co-maintain auto-regulation

Many studies indicate that PTF1 complex involves in the formation of acinar cells during pancreatic development. The protein Ptf1a cooperates with one of the common E-proteins to bind RBPJ/RBPJL to produce the PTF1 complex PTF1-J/PTF1-L, which can form auto-regulatory loops to maintain the expression of Ptf1a. The details of the regulation are mentioned in Section 3.1.2, and the auto-regulatory loop of Ptf1a via PTF1 complex is shown in Figure 8(A). Meanwhile, some studies suggest that Mist1 involves the formation of auto-regulatory loop related to Ptf1a in acinar cells ^[37,44]. Ptf1a drives transcription of Mist1, which induces expression of the connexin, and it in turn promotes the formation of Ptf1a, thus forming an auto-regulatory loop via Mist1 which maintains the acinar cell phenotype ^[44], as shown in Figure 8(B). In addition, PTF1 complex induces transcription of Mist1, which results in the maintenance of Ptf1a auto-regulation ^[37]. Moreover, the loss of functional Mist1 results in progressive acinar damage and the acquisition of certain ductal properties ^[45]. Perturbations to the complex PTF1-L affect the phenotype of the mature acinar cells, i.e., low differentiation of acinar tissue ^[36]. These results indicate that both PTF1 complex and Mist1 play important roles in maintaining formation of the auto-regulatory loop related to Ptf1a, and perturbations to one of them affect the formation of acinar cells, as shown in Figure 8(C). Therefore, the auto-regulatory loop involving Ptf1a in Figure 8(C) can explain the phenomenon that inhibition of PTF1 complex or Mist1 expression leads to damaged acinar cell formation because inhibition of PTF1 complex or Mist1 can break the Ptf1a auto-regulatory loop.

3.2.   The regulatory function of Notch signaling in pancreatic cell differentiation

The communication among cells are necessary for homeostasis in multicellular organisms. In fact, intercellular Notch signaling plays a key role in the decision between endocrine and exocrine cells. The Notch intracellular domain (NICD) is generated by intercellular interaction and then inhibits the expression of Ngn3 indirectly ^[20]. Studies suggest that intercellular communication via Notch signaling, which can be cis-inhibited by Dll1 within the same cell, and trans-activated by Dll1 in adjacent cells ^[24], plays a vital role in the pancreatic cell destiny decisions ^[4,18,46]. The pro-endocrine protein Ngn3 activates the expression of Dll1 ^[18]. Notch binds to Dll1 of adjacent cells, leading to the release of NICD. NICD activates expression of Hes1 which represses Ngn3 ^[44]. NICD regulates Ngn3 and Ptf1a expression, i.e., down-regulation of Notch results in increased expression of Ngn3 and decreased expression of Ptf1a ^[18]. These results indicate that Notch signaling affects fate decision of endocrine/exocrine progenitor cells in the first stage of differentiation. Indeed, Notch signaling is induced in exocrine progenitors and down-regulated in endocrine progenitors ^[47]. However, Notch signaling suppresses liveness of the Ptf1 transcriptional complex (PTF1) ^[47], and the PTF1 complex is necessary for exocrine cell differentiation, e.g., acinar cells. Therefore, Notch signaling also plays a vital role in the transformation of exocrine progenitor into exocrine (acinar) cells in the second stage of differentiation.

3.2.1.   Regulation of Notch signal in the first stage of differentiation

For the sake of simulation, we consider the interaction between only two cells, i.e., $N = 2$ in Eq (2). In addition, we set $a_1 = a_2 = a_{NP}$ with $a_{NP}$ as the control parameter. All values of other parameters are given in Table 4. The change of Notch signaling is realized through change of the basal production rate of Notch, i.e., high level at $\beta_{N, i} = 20$, $\beta_{D, i} = 20$, and low level at $\beta_{N, i} = 10$, $\beta_{D, i} = 20$. Now, we explore the regulation of Notch signaling on endocrine/exocrine progenitor cell fate decisions from the perspective of bifurcation analysis. At the low level, cells stay at the pluripotent progenitor state with low Ngn3 $(x_{2, i})$ and Ptf1a $(x_{3, i})$ in cell $i$ at small $a_{NP}$, i.e., $a_{NP} < a_{{NP}_{c_1}}$. As $a_{NP}$ increases, pluripotent progenitor cells differentiate finally into endocrine progenitor cells with high Ngn3 $(x_{2, i})$ and low Ptf1a $(x_{3, i})$ in cell $i$ at $a_{NP} > a_{{NP}_{c_1}}$, as show in Figure 9. While at the high level, when $a_{NP} < a_{{NP}_{c_2}}$, cells also stay at the pluripotent progenitor state. As $a_{NP}$ increases, pluripotent progenitor cells differentiate finally into exocrine progenitor cells with high Ptf1a $(x_{3, i})$ and low Ngn3 $(x_{2, i})$ in cell $i$ at $a_{NP} > a_{{NP}_{c_2}}$, as show in Figure 10.

These results show that down-regulation of Notch signaling may induce committed endocrine progenitors. While relatively strong Notch signaling leads to cell differentiation into the exocrine progenitors state, as shown in Figures 9 and 10. In fact, in the Notch signaling, NICD interacts with the DNA-binding protein RBP-Jk to induce expression of the bHLH $Hes$ genes, hereby inhibiting the expression of downstream target genes involved $Ngn$ genes. The results in ^[18] showed that absence of Dll1 or RBP-Jk (impaired Notch signalling) can result in increased expression of the pro-endocrine gene $Ngn3$, thus promoting the endocrine destiny. However, at the normal level of Notch signaling, cells express Hes1 and p48 and select the direction of exocrine destiny. In addition, nuclear Hes1 protein is also present in the most cells expressing Ptf1-p48, suggesting that Notch signaling is activated in an exocrine progenitor pool ^[47].

3.2.2.   Regulation of Notch signal in the second differentiation stage

The roles of Notch signaling on endocrine/exocrine differentiation in the first stage is mentioned in Section {3.2.1}. In fact, Notch signaling also plays a vital role in decision fate of exocrine cells during the second differentiation stage, e.g., absence of Notch function is related to accelerated development of acinar cells during developing pancreas ^[47,49]. In addition, Notch signaling suppresses liveness of the Ptf1 transcriptional complex (PTF1), resulting in cells remaining in the exocrine progenitor state ^[47].

We here just propose a network shown in Figure 11(A) by combining the simplified network in Figure 2 with the auto-regulatory loop of Ptf1a in shown Figure 8(C). Further modeling and analysis of Figure 8(C) will be performed in future. In addition, it is assumed that NICD inhibits PTF1 ^[50]. In fact, it can be seen that inhibition of PTF1 activity by NICD can hinder the formation of auto-regulatory loop of Ptf1a, according to the analysis of network in Figure 11(A). It is mentioned in Section {3.1.2} that the auto-regulation of Ptf1a is necessary for the formation of acinar cells, so the up-regulation of Notch signaling (NICD) inhibits the differentiation of acinar cells, thus maintaining the state of exocrine progenitor cells in the second differentiation stage. Although Notch signaling inhibits the activity of PTF1, the details of how Notch signaling regulates PTF1 is still not clear. It is possible that Notch signaling inhibits the expression of E-cadherin, leading to breaking of the auto-regulation loop of Ptf1a and inhibition of PTF1, as show in Figure 11(B). The loss of Notch function during the development of zebrfish pancreas is related to the accelerated development of exocrine pancreas, which may be related to a mechanism that Notch inhibits the activity of Ptf1 complex (PTF1), independent of the change of Ptf1 component protein levels, and thus affects the development of exocrine pancreas ^[47]. According to the model in Figure 11(B), loss of Notch can reduce the amount of NICD, thus promoting the expression of PTF1 to realize the fate decision of exocrine cells. Moreover, Ptf1a-expressing cells participate in lateral stabilization in pancreatic cell development, such conditional activation is, in principle, consistent with both cadherin/beta-catenin signaling. When cadherin-mediated cell-cell adhesion is disrupted, lateral stabilization is lost, leading to acinar-to-islet cell conversion ^[22]. In fact, it can be inferred from the model in Figure 11(B) that if E-cad expression is inhibited in acinar cells, Ptf1a is difficult to maintain, resulting in acinar transform into islet cells.

3.2.3.   Lateral inhibition regulates the fate decision of adjacent pancreatic cells

Lateral inhibition is a process by which cell communication is realized so as to drive neighboring cells to take different destinies. Notch signaling is a typical mediator of lateral inhibition ^[23]. In fact, a small difference in the generation rates of Notch and Delta results in a much larger difference in Notch signaling activity, thereby promoting the formation of lateral inhibition patterning ^[48]. In addition, lateral inhibition via Notch signaling plays important roles in maintaining endocrine cell scattering distribution ^[18,44]. Although we know that lateral inhibition is related to the regulation of pancreatic cell destiny decisions, little is known about how it is regulated. Now, we explore the mechanism of lateral inhibition in pancreatic fate decisions by bifurcation analysis of mathematical model Eq (2).

For convenience, we consider the interaction via Notch signaling between only two cells, i.e., $N$ = $2$ in Eq (2). In addition, we set $a_1$ = $a_2$ = $a_{NP}$ with $a_{NP}$ as the control parameter. Especially, the following parameter values are used as standard values $\beta_{N, 1}$ = $5$, $\beta_{D, 1}$ = $20$, $\beta_{N, 2}$ = $20$, $\beta_{D, 2}$ = $20$. Other parameter values are given in Table 4. The difference in generation rate of Notch ($\beta_{N, i}$) leads to the difference in Notch signaling activity between the two cells. When the difference is large enough, two adjacent cells may adopt different cell fates. In other words, one pancreatic progenitor cell with low Ngn3 $(x_{2, 1})$ and low Ptf1a $(x_{3, 1})$ at $a_{NP} < a_{{NP}_{s_1}}$ can differentiate into either a endocrine progenitor cell with high Ngn3 $(x_{2, 1})$ and low Ptf1a $(x_{3, 1})$ at $a_{NP} > a_{{NP}_{s_1}}$, as shown in Figure 12(A), (B), while another pancreatic progenitor cell with low Ngn3 $(x_{2, 1})$ and low Ptf1a $(x_{3, 1})$ at $a_{NP} < a_{{NP}_{s_2}}$ may differentiate into a exocrine progenitor cell with low Ngn3 $(x_{2, 1})$ and high Ptf1a $(x_{3, 1})$ at $a_{NP} > a_{{NP}_{s_2}}$, as shown in Figure 12(C), (D). These results shows that lateral inhibition regulates the fate selection of exocrine/endocrine cells by controlling the first differentiation stage. The bifurcation diagram of Figure 12 analysis show that our mathematical model can explain the model of lateral inhibition from Notch-Delta1 in Cell1 and Cell2 fate of Figure 13.

Since our current model is built on intertwined dynamics with trans-activation and cis-inhibition involving the Notch and Delta proteins ^[24,48], therefore, except for the pancreas cell fate decisions involved in Notch-Delta of lateral inhibition mechanism, there are many other cell fate decisions involved this mechanism, including angiogenesis ^[53], spinal cord patterning in zebrafish ^[54], and development of neuroblast cells in early neurogenesis ^[55], can be explained by our current model. In addition, we do not consider the function of cis-inhibition in the regulation of pancreatic cell fate decision in our study. In fact, Delta1-mediated cis-inhibition is necessary to control cell fate selection ^[25]. From the analysis of our regulatory network in Figure 2 in the first stage of differentiation, it may be deduced that when Ngn3 in progenitor cells stimulates the expression of Delta1, increased Delta1 will cause Notch to be cis-inhibited, resulting in less Notch signaling, thus maintaining the expression of Ngn3 to reach the threshold of endocrine progenitor commitment.

Remark 2. Different levels of NICD in the two cell system induce two different bifurcation points. When two cells are coupled, each cell has two saddle-node bifurcation points $a_{{NP}_{s_1}}$ and $a_{{NP}_{s_2}}$, as shown in Figure 12. Actually, the bifurcation point $a_{{NP}_{s_1}}$ is induced by the first cell, while the bifurcation point $a_{{NP}_{s_2}}$ is induced by the second cell. As a whole, the system has two bifurcation points. Due to the differences between two cells, they adopt different cell fates under two different bifurcation points. More exactly, the first cell transforms into an endocrine progenitor cell, while the second cell becomes an exocrine progenitor cell.

Remark 3. In fact, the Ngn3 in pancreatic endocrine progenitor cells controls the expression of the Notch ligand Delta1, which activates the expression of Notch target genes such as Hes1 and thereby represses endocrine differentiation in neighboring cells (Notch trans-activation) by lateral inhibition ^[49,52]. Moreover, the Ptf1a via activation of Delta1 stimulates multipotent pancreatic progenitor cells proliferation and contributes to Hes1 activation, and thereby may indirectly contribute to maintaining high Ptf1a protein levels ^[52].

4.   Conclusions

Our work enriches the research on fate selection of pancreatic cells and help people to understand the mechanism of pancreatic cell development from multiple perspectives. So far, there are few researches related to pancreatic cell development by using mathematical models. We construct related regulatory networks at the levels of single cell and multiple cells, and further build mathematical models and analyze them by bifurcation analysis, which helps us understand the mechanism of pancreas cell fate decisions more intuitively. Meanwhile, we divide the process of exocrine/endocrine cell differentiation into two stages and reveal the regulatory mechanisms of auto-regulation loop of Ptf1a and acinar-to-islet cell conversion. It is known that auto-regulation of Ptf1a is essential for the formation and remain of mature acinar cells, but the specific regulation mechanism of Ptf1a auto-regulation is still not clear. So, we provide related information about the auto-regulation of Ptf1a by extensive literature studies ^[30,37,44]. It can help us understand the phenomenon of damaged acinus cell formation caused by inhibited expression of PTF1 complex or Mist1 expression. Although the mechanism of pancreatic cell fate regulation was explored in ^[15], but the detailed regulatory mechanism of exocrine progenitor cell differentiation has not been considered. Here, we further consider detailed regulatory mechanism of exocrine cell differentiation into acinar cells in our study. The ^[30] revealed that auto-regulation loops of Ptf1a and Rbpjl may play an important role in maintaining a stable phenotype of pancreatic acinar cells, but the mechanism of the auto-regulation loop in pancreatic cell fate decision has not been analyzed from the perspective of mathematical models. Compared to the literature ^[30], we qualitatively analyze the mechanism of the auto-regulation loops in the formation of acinar cells from the point of view of mathematical model.

Notch signaling takes an important part in the process of exocrine and endocrine development mainly at two stages of differentiation. At the first differentiation stage, high and low level of Notch signaling forces pancreas progenitor cells to differentiate into the exocrine progenitor and endocrine progenitor lineages, respectively. At the second differentiation stage, Notch signaling inhibits the formation of acinar cells so that cells remain in the state of exocrine progenitor lineage, which is induced because NICD inhibits the expression of the PTF1 complex and breaks the auto-regulation of Ptf1a. The mechanism by which NICD down-regulates PTF1 complex expression is still poorly studied. In order to further understand the regulatory mechanism, we propose a possible regulatory network in Figure 11(B) according to relevant studies ^[44,50,51]. The regulatory role of Notch signaling in pancreatic cell differentiation was explored in ^[18], and lateral inhibition mechanism based on Notch signaling was mentioned in ^[49]. However, compared with references ^[18,49], we enrich and develop mathematical model to understand the role of Notch signaling and lateral inhibition in the regulation of pancreatic cell fate.

Acknowledgments

This research is supported by the National Natural Science Foundation of China (Grant No. 11971297).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix.  

Bifurcation theory analysis of Eqs (1) and (2)

Firstly, we analyze these bifurcation conditions in Figure 3(A) based on nonlinear dynamic Eq (1). At $n = 2$ and $n_1 = 6$, the steady state $(\varphi_1, \varphi_2, \varphi_3, \varphi_4, \varphi_5)$ of Eq (1) satisfies

In addition, the Jacobian matrix of Eq (1) is

Where

The characteristic equation of the Jacobian matrix (4) can be obtained $|\lambda E-J| = 0$, bifurcation occurs when at least one eigenvalue is equal to 0, it has

The bifurcation point of Eq (1) can be obtained by combining Eqs (3) and (5). For Figure 3(A), the $a_1 = a_2 = a_{NP}$ as the control parameter, the other parameters are predetermined values in Table 3. The steady state of Eq (1) under this parameter condition can be described as

It has $(\varphi_4, \varphi_5) = (0, 0)$ or $(\frac{5\pm\sqrt{23}}{10})$ by analyzing Eq (6). And considering that Isl1 ($x_4$) and Rbpji ($x_5$) are not activated in the case of Figure 3(A), i.e., its expression level is zero, $(\varphi_4, \varphi_5) = (0, 0)$. We know the value of $\varphi_4$ and $\varphi_5$, so, it has $\varphi_1 = 1$. The Eq (7) can be obtained by simplifying Eq (6) under $\varphi_1 = 1, \varphi_4 = 0, \varphi_5 = 0$.

Moreover, we substitute the parametric values Figure 3(A) and partial equilibrium ($\varphi_1 = 1, \varphi_4 = 0, \varphi_5 = 0$) of Eq (6) into Eq (5). It has Eq (8), as following:

We can obtain the bifurcation condition in Figure 3(A) and find the bifurcation point by combining Eqs (7) and (8). Similar to the above analysis, the conditions for their bifurcation can be analyzed by substituting the parameter values in Figures 4–6 into Eqs (3) and (5), respectively.

Then, we analyze these bifurcation conditions based on nonlinear dynamic Eq (2) and assume the two cells $(i = 1, 2)$ that communicate are the same cells, $x_{k, 1} = x_{k, 2} = \varphi_k (k = 1, 2, ..., 8)$ and the parameters $\beta_{D, 1} = \beta_{D, 2} = \beta_D, \beta_{N, 1} = \beta_{N, 2} = \beta_N$, it can reduce the dimension of the system to facilitate analysis. At $n = 2$ and $n_1 = 6$, the steady state $(\varphi_1, \varphi_2, \varphi_3, \varphi_4, \varphi_5, \varphi_6, \varphi_7, \varphi_8)$ of Eq (2) satisfies

Naturally, we can get the Eq (9) Jacobian matrix

Where, the formula for $\phi_1, ..., \phi_{12}$ are the same as Eq (4). $\phi_{13} = -\frac{2a_1 b \varphi_1^2 \varphi_8}{K1^2+\varphi_1^2+b_3 \varphi_3^2+b \varphi_8^2}$,

$\phi_{14} = \frac{2a_8 K_9^2 \varphi_2}{(K_9^2+\varphi_2^2)^2}$, $\phi_{15} = -\varphi_7(\frac{1}{k_c}+\frac{1}{k_t})$, $\phi_{16} = -\varphi_6(\frac{1}{k_c}+\frac{1}{k_t})$, $\phi_{17} = \frac{\varphi_7}{k_t}$, $\phi_{18} = \frac{\varphi_k}{k_t}$.

The characteristic equation of the Jacobian matrix (10) can be obtained $|\lambda E-A| = 0$, bifurcation occurs when at least one eigenvalue is equal to zero, it has

By substituting the parameter values in Figures 9 and 10 into Eqs (9) and (11), the bifurcation points in Figures 9 and 10 can be obtained, respectively. In addition, Eq (2) under the condition of Figure 12 will be more complex, because part of the parameters of the two cells are different, so it is not possible to simplify as Eq (2) under the condition of Figures 9 and 10. Although it is a 16-dimensional system, the theoretical analysis method is the same as above.

Table 3: Here, $a_1$ = $a_2$ = $a_{NP}$ as control parameter and $a_3$ = $a_4$ = $a_6$ = $0$, $a_7$ = $50$, $b_3$ = $400$ in Figure 3; $a_1$ = $a_2$ = $2$, $b_3$ = $400$, $a_3$ = $a_7$ = $50$ and $a_4(a_6 = 5)$, $a_6(a_4 = 2)$ as control parameter in Figures 4 and 5; $a_3$ = $a_7$ = $a_P$ as control parameter and $a_1$ = $a_2$ = $2$, $a_4$ = $2$, $a_6$ = $5$, $b_3$ = $400$ in Figure 6; $a_3$ = $a_7$ = $a_P$ as control parameter and $a_1$ = $a_2$ = $2$, $a_4$ = $2$, $a_6$ = $5$, $b_3$ = $50$ in Figure 7. Some parameter values in the table are extracted from reference ^[22], i.e., $\alpha$ = $1$, $a_5$ = $50$, $b_1$ = $100$, $b_2$ = $100$, $d_i$ = $1$ $(i = 1, 2, 3, 4, 5)$, $K_i$ = $1$ $(i = 1, 2, 3, 5)$ and $n = 2$. However, the remaining parameter values in the table are obtained based on the parameter values in which bifurcation behavior can occur in Eq (1).

Table 4: Some parameter values in this table were obtained by referring Table 3 and supplementary information of reference ^[48], i.e., $d_8$ = $1$, $k_c$ = $1$, $\beta_{N, i}$ = $\{5, 10, 20\}$ and $\beta_{D, i}$ = $20$. The remaining parameter values in the table are obtained based on the parameter values in which bifurcation behavior can occur in Eq (2).