Mapping known IL-1 regulations to a mathematical modelIn order to minimize the number of unknown parameters and variables, we aimed to construct a simple model, which captures the qualitative behavior of the biological system. In our earlier theoretical studies, we have shown that cytokine auto-stimulation through the NF-κB pathway belongs to a class of phenomena known as excitable media27,28,29. Here we modify an earlier model27,29 to describe the NF-κB mediated inflammatory response in pancreatic islets. It has been shown that, in addition to macrophages, β-cells contribute to IL-1β accumulation in the islets30. In the current model, we do not differentiate between β-cells and macrophages, and refer to them as IL-1β Secreting Cells (ISC). We assume that other endocrine cell types (α, δ and γ) do not contribute to cytokine secretion.A two-dimensional model was constructed in which cells were arranged on a hexagonal grid, and their positions were fixed. All cells inside a circular area in the middle of the hexagonal grid were modeled as ISCs representing a pancreatic islet (Fig. 2B). Cells outside the islet represented surrounding tissue in which IL-1β was allowed to diffuse and degrade. IL-1β was assumed to diffuse freely between cells; based on a molecular weight of 17 kDa, the diffusion coefficient of freely diffusing IL-1β is estimated to be D = 20 µm2/s31. While we do not explicitly model small blood vessels penetrating into islets, we account for the fact that β-cells are typically in close vicinity to the blood, by using plasma clearance rates of IL-1β as an estimate of the effective IL-1β half-life in the islets \({\boldsymbol{(}}{{\boldsymbol{\tau }}}_{{\boldsymbol{I}}}\,{\boldsymbol{=}}\,{\mathbf{0}}{\mathbf{.}}{\mathbf{17}}\,{\boldsymbol{h}}{\boldsymbol{)}}\)32. The choice of these and other parameters are summarized in Table 1, “Methods”.Fig. 2: Transient and persistent IL-1 modes emerge in single islets.A Snapshots of simulated IL-1 dynamics in an islet in response to different combinations of high and low S1 (Signal 1) and S2 (Signal 2). All cells within the red circle are modeled as ISCs (intra-islet cells), and outside of the circle as non-secreting cells (surrounding tissue cells). Only the S1 sources cells (red dots) have non-zero values of S1 between 0 and 10 h (B, top panel). B Time courses of IL-1β averaged over all cells in the islet corresponding to the islets in (A). The gray area indicates the time interval in which the S1 sources are “tuned on”. Under all conditions, the islet displays an initial strong secretion of IL-1β peaking after ~2 h. At low S1 and S2 (top panel), the islet responds transiently and will return to the resting state (Supplementary Movie 1). If either S2 or S1 is high (middle and bottom panel), the islet enters a persistent (“locked”) state with sustained IL-1β expression (Supplementary Movies 2 and 3). C, D Probability of locking at different levels of S1 and S2. The random position of S1 sources (indicated by red dots in A) introduces stochasticity in the islet fates (locking) for some values of parameters S1, ns and S2 (see also Fig. 5D). Here, we show the results of 50 simulations with equal parameter values, but with different random positions of the S1 sources. The fraction of persistent islets increases with the total S1 dose (number of sources (ns) times the S1) (C) and with increasing S2 (D). Gray points and error bars show the average number of persistent islets, gray line is a sigmoidal fit to the fraction of persistent islets, used to find EC50 (see “Methods”, Supplementary Note 4 and Supplementary Fig. S5). Overlaid in blue and yellow are the bee-swarm plots, where each dot represents five simulations (all yellow dots have a value of 1 and all blue dots have value of 0). The bee-swarm plots visualize that islets are in either transient (blue dots) or persistent (yellow dots) states, as indicated by the insert snapshots to the far right. For both parameter scans (C and D), there is a region of co-existence between persistent islets and transient islets that return to the resting state. The state of individual islets depends on the positions of S1 sources. Here each islet has a radius of 7 ISCs. Results for larger sizes are shown in Fig. 5.The individual ISCs were modeled as entities that were able to both sense surrounding levels of IL-1β and to secrete IL-1β into the extracellular space. We reduced the detailed regulatory network of IL-1β, shown in Fig. 1A, to an effective regulatory network (Fig. 1B). The model includes four rescaled dynamic variables: I, I’, N and R (see details about rescaling in “Methods”). I represents the concentration of IL-1β, I’ is the concentration of pro-IL-1β, N represents the concentration of nuclear NF-κB (i.e., active in transcription), and R represents the combined effect of regulating proteins that inhibit active NF-κB (e.g., IκBα, A20, cylindromatosis (CYLD), etc.). While the IκBα negative feedback is significantly faster, compared to A20 and CYLD, combining them together in one slow feedback produces qualitatively the same results as long as one does not aim to reproduce the lower amplitude secondary oscillations27.The model for an ISC is therefore a system of ordinary differential equations together with a partial differential equation for IL-1β:$$\frac{{\rm{d}}N}{{\rm{dt}}}={k}_{{IN}}\frac{{I}^{2}}{{I}^{2}+1}\left(1-N\right)-\,{k}_{{RN}}R$$
(1)
$$\frac{{\rm{d}}R}{{\rm{dt}}}={k}_{{NR}}N-\,R/{\tau }_{R}$$
(2)
$$\frac{{\rm{d}}I^{\prime} }{{\rm{dt}}}={pN}-\left({S}_{2}+\frac{1}{{\tau }_{{I}^{{\prime} }}}\right)I^{\prime}$$
(3)
$$\frac{{\rm{d}}I}{{\rm{dt}}}={{S}_{1}\delta (\left\{{x}_{i}\right\})+S}_{2}{I}^{{\prime} }-\frac{I}{{\tau }_{I}}+D{\nabla }^{2}I$$
(4)
The two-dimensional model of a pancreatic islet is a reaction-diffusion model on a hexagonal grid. While the variable I is allowed to diffuse between cells, the variables N, R and I’ are local intracellular variables. Only the cells representing ISCs in the pancreatic islet are able to synthesize de novo IL-1. This is modeled by setting the parameter p to a positive value for ISCs cells and to zero for all non-ISCs.The first two equations governing dynamics of NF-kB activation and inhibition (Fig. 1D) focus on capturing the characteristic initial peak in NF-kB response and were first introduced in our earlier work27,29.The parameters kIN, kRN, kNR and τR have been fixed by fitting the model to experimental data for TNF-induced NF-κB dynamics in single cells33, such that the timing of the initial peak in NF-κB activation matched experimentally observed values27,29,33 (see “Methods” for details of fitting).While the data on IL-1β induced NF-κB responses in single cells is not available; population-averaged data supports transient activation of NF-κB with a peak around 2 h34.The last two equations and parameters p, S1, S2, τI, τI’, and D are related to the dynamics (production, degradation and diffusion) of IL-1β. Our investigation focuses on the parameters S1 and S2, related to the levels of initiating inflammatory (Signal 1) and metabolic cues (Signal 2), and how these parameters affect the dynamics of the model. Further details on model parameters and approximations are presented in detail in Table 1 in the “Methods” section.Although IL-1β could be produced solely by islet-invading immune cells such as macrophages, it has been hypothesized that the initial increase of IL-1β is mainly produced by stressed β-cells. This increase may subsequently attract immune cells, adding further to the increasing levels of IL-1β14. In our model, we do not specify the originators of the initial increase of IL-1β, but simulate Signal 1 by adding an additional production of IL-1β (the term \({S}_{1}\delta (\left\{{x}_{i}\right\}\) in Eq. (4)) to a number of randomly chosen ISCs within the islet. These cells thus become sources of Signal 1. As these cells are chosen at random, the positioning of these Signal 1 sources will be also random for each simulation. In Fig. 2A and throughout, the position of Signal 1 sources are indicated by red dots. The sources of Signal 1 act as initiators of the model, allowing us to study how the ICSs respond and regulate IL-1β when exposed to an external stimulus. While all the ISCs are able to produce and secrete IL-1β, only Signal 1 sources have non-zero value of the S1 parameter. Overall Signal 1 in the system is thus characterized by both the strength of basal IL-1β production, S1, and the number of sources, ns (total Signal 1 dose is equal to \({S}_{1}\times {ns}\)). We stress that the initial sources described by the S1 parameter reflect a cumulative effect of other NF-κB activating agents such as TNF, TLR, PAMPs or other cytokines, but for the purpose of our model, we describe these as IL-1β sources in order to minimize the number of variables. While we choose to model S1 as a few discrete sources, we have tested that the model results hold if, instead, S1 is equally distributed across all cells. This configuration is, in particular, relevant for cases where Signal 1 represents low-grade inflammation stemming from sources outside the islets, with low levels of pro-inflammatory cytokines circulating in the blood35.The severity of T2D is positively correlated with the levels of glucose, the amount of free fatty lipids concentration in the blood and DAMPs in the extracellular space. These factors all contribute to the activity of the NLRP3 inflammasome and, in turn, caspase-1, see Fig. 1A. In order to construct a simple model, we merge the effects which contribute to an increased caspase-1 activity into the S2 parameter (see Fig. 1C). In our simulations S2 is set to be the same for all ISC and does not change in time. But since caspase-1 actively cleaves pro-IL-1β into mature IL-1β, the closest biological interpretation of S2 is the activity of caspase-1. S2, in combination with the production rate of pro-IL-1β, p, constitutes the strength of the positive feedback by which β-cells amplify the local IL-1β concentration. This becomes a crucial parameter in the qualitative behavior of the simulated islets.The ISCs synthesize IL-1β and secrete it into the extracellular space, where IL-1β acts in an autocrine and paracrine manner. In each simulation, islets are initiated in a steady state of low activity with no stimulation, i.e., NF-κB is inactive and S1 equals zero for all ISCs. At time zero the islets are stimulated by a number of small inflammation sources—technically S1 is set equal to a non-zero value for a discrete set of spatial positions. The sources are only non-zero for a finite time, after which the S1 parameter is set equal to zero, and we monitor the islet response. While the S2 parameter is non-zero at all times and for all ISCs (cells within the islet, indicated by a red circle in Fig. 2A), the S1 parameter is only non-zero for a discrete set of islet cells (indicated by red dots in Fig. 2A).The model predicts transient and persistent production of IL-1β in single islets, in correspondence with the IL-1β dual effectThe bimodal effects of IL-1 on β-cell function and viability are associated not only with the concentration of IL-1 but also with the duration of exposure. Transient exposure leads to protective effects, while prolonged exposure promotes cell death14,36,37. As pancreatic islets have been shown to mount a significant local IL-1 response, among the endocrine islet cells and infiltrating macrophages, we wanted to investigate if the two distinct temporal profiles of IL-1 may emerge at the single islet level. To investigate this, we monitored how levels and duration of the IL-1 production depend on inflammatory Signal 1 (S1) and nutritional Signal 2 (S2) cues within an in silico islet.We found that whereas low levels of S1 and S2 result in transient upregulation of IL-1β, high levels of either S1 or S2 (or both) will result in sustained high levels of IL-1β (Fig. 2A, B). When an islet with a low S2 is stimulated by a small number of inflammatory sources (low ns), it responds with a pulsating NF-κB activity and a similar pulsating amplification of IL-1β (Fig. 2A). Once the sources are removed, the islet stops pulsing and returns to the resting state (see Supplementary Movies 1 and 4). Hence, the islet exhibits a transient response to the transient inflammatory cues. If a similar islet is exposed to either higher Signal 1 (more initial sources of inflammation) or to higher Signal 2 (e.g., higher caspase-1 activity), the cells will collectively amplify the IL-1β concentration and transition into a ”locked” state of constantly elevated IL-1β (see Supplementary Movies 2, 3, 5 and 6 and Supplementary Note 1). In this state, the islet will sustain the high levels of IL-1β, and will not settle back to the resting state—even when the initial IL-1β sources are removed. The locked state simulates a state of chronic inflammation, where prolonged exposure to high levels of IL-1β has been reported to have deleterious effects leading to impaired insulin secretion and β-cell death14,15,16,17. Similarly to our earlier work,27,29, dynamical systems analysis showed that the locked state is a consequence of an overactive positive feedback in ISCs, which prevents the fast variables (I, I’ and N) from undergoing a saddle node bifurcation and return to low cytokine levels. These results are consistent with the dual role of IL-1β reported, where it can both facilitate cell survival, the transient response, or enter a state of chronic inflammation, the locked state14,17. Notably, the islets display at least one NF-κB/IL-1β pulse in response to a stimulus, and therefore, even the islets that eventually lock also display an initial transient IL-1β response (Fig. 2A, B), suggesting that the chronic inflammation state will always follow after the initial acute phase with transient IL-1β.Having established two distinct modes of IL-1β response, we investigated if the transition between the two modes is gradual, with, e.g., islets gradually increasing production levels of IL-1β with increasing S1,2 or sudden, with a discontinuous jump from low to high IL-1β.We found that the transition is sharp (Fig. 2C, D and confirmed by the larger parameter scan shown in Fig. 3E–H): a small increase in either Signal 1 or Signal 2 (or both) can have a drastic effect on the probability of locking. Notably, close to the transition, locked and non-locked islets can coexist: despite being exposed to the same values of S1, ns and S2, some islets lock into a state with sustained level of IL-1β, whereas others show transient behavior. In this “co-existence regime,” the fate of an individual islet is determined by the exact geometrical positioning (and timing) of the S1 sources. This finding is in line with the current knowledge that in the same pancreas there can be co-existing subpopulations of inflamed and non-inflamed islets15 (measured by number of infiltrating macrophages).Fig. 3: Transition from IL-1 protective to deleterious effects is accelerated under higher S2.A–D Experimental data for the time course of the Immunoreactive Insulin Release (IRIR) normalized with respect to IRIR in absence of IL-1β based on data from ref. 37. Color codes for the fold change in IRIR in presence of the IL-1 relative to the IRIR in absence of IL-1. The doses of IL-1 (corresponding to S1) are specified on the x-axis, and glucose (corresponding to Signal 2, S2) levels on the y-axis. A Islets were exposed to the specified glucose and IL-1 levels for 6 h, and the color indicates measured normalized IRIR over the period 0–6 h. B Similar to (A), but cumulative IRIR was measured for the period of 6–24 h. C Similar to (A), but cumulative IRIR was measured for the period of 24–48 h. D Similar to (A), but cumulative IRIR was measured for the period of 48 h–6 days. I IRIR in absence of IL-1β data from ref. 37. The IRIR was measured in vitro in islets exposed to different concentrations of glucose and IL-1β. The islets were exposed to increasing durations of IL-1β from (A) to (D). E–H Simulation results showing the fraction of simulated islets responding transiently at a given value of Signal 1 dose, S1 x NS (x-axis)and S2 (y-axis). The fractions are color-coded from yellow corresponding to 0, and thus no islets in transient state (i.e., all are locked) to blue corresponding to 1 (all islets responding transiently). The fraction for each circled point is calculated based on 20 simulations of individual islets (similar to Fig. 2A) of a radius of 6.7 cells with randomly placed S1 sources within the islet. Similar to the experiments, the probability of transitioning to the “locked state” increases as function of external S1 stimulus, inflammasome activation (S2) and the duration of external S1 stimulus. Each external source has an S1-value of 25 h−1 and p = 1400 h−1, NS is a number of sources.The model recapitulates the fact that glucose narrows the time window in which IL-1 has beneficial effectsHaving established the emergence of transient and persistent modes of IL-1 production within the single islet, we wanted to further test if our model is consistent with the experimental observation that the time window of the beneficial effect from IL-1 is shorter at higher glucose levels (corresponding to higher S2)37.As mentioned in the “Introduction”, the beneficial effect of transient IL-1 exposure is manifested by increased insulin release. However, long exposure to high levels of IL-1 leads to cell death and decreased insulin production within the islet. Therefore, the ability of the islet to produce insulin can be used as a proxy of its state. For simplicity, we assume that the in silico islet locked in a persistently elevated IL-1 state is representative of the real islets where IL-1 exposure has deleterious effects, i.e., those with decreased insulin release.Palmer et al.37 measured the immunoreactive insulin (IRI) release in islets that had been exposed to different glucose concentrations (3.3, 5.5 and 11 mM) for different durations of time (6, 24, and 48 h and 6 days of culture). As insulin release increases with increasing glucose concentrations in the absence of IL-1β (Fig. 3I), to extract the effects of IL-1, one has to normalize for these baseline IRI releases. To test the complex dual effects of IL-1β on β-cell function different concentrations of IL-1β (50, 500, and 2000 ng/L) was added to the culture, and IRI release was measured. To focus on the effects of IL-1β, we have plotted the values reported in ref. 37 in Fig. 3A–D, where we normalized IRI release at each glucose and IL-1β level by the corresponding values in the absence of IL-1β (Fig. 3I). Qualitatively, the effects of IL-1 are glucose-independent at brief (0–6 h, Fig. 3A) and long (2–6 days, Fig. 3D) exposures. At all glucose levels, brief exposure to IL-1 increases, and long exposure to IL-1 inhibits insulin release. The switch between the two modes, which happens at intermediate IL-1 durations, is however glucose-dependent, occurring already after 6 h at the highest glucose concentration (Fig. 3B).For our model to be consistent with these results, brief exposure to S1 should result in a majority of islets in transient mode, which is indeed the case (Fig. 3E, F). On the other hand, experiments showed that long exposures to IL-1β (48 h–6 days) lead to decreased IRI release when compared to the reference values (Fig. 3D). This is also consistent with our model outcomes for the situation where the S1 sources are turned on for an increased period of time (1 h, Fig. 3H). In this case, the model predicts increased probabilities of islets locking into the persistent mode, which we interpret as being related to impaired insulin response. We note that the exact timing of the model does not match the longer time period investigated in the experiments.At the intermediate durations (6–24 and 24–48 h), experiments show that switching from the transient to the persistent mode should happen first where S2 is highest (Fig. 3B, C). We find that this is indeed the case (Fig. 3E–H). Notably, the islets exposed to IL-1β for a longer time display an increasingly impaired response to glucose stimulation. This is consistent with our simulations, where the parameter range of ns, S1 and S2 that leads to locked islets increases with longer exposure to S1.To our knowledge, it is not known why the switching from insulin-enhancing to insulin-inhibiting modes happens earlier at higher glucose concentrations. While it has been shown that glucose potentiates IL-1β-induced nitric oxide production—one mediator of the IL-1β cytotoxic effect— the time dependence has not been addressed38. Our results suggest that the time dependence may come as a result of the modulation of the positive feedback loop by S2, in that stronger positive feedback accelerates the transition to “locked” state and thus to insulin-inhibiting mode of IL-1β. When we compare simulation and experimental results, the region of parameters where we find locked islets increases much faster in simulations than the corresponding region for the islets with impaired Immunoreactive Insulin Release (IRI) release in the experiments, Fig. 3. In the experimental settings, inflammasome activation (Signal 2) is likely to increase gradually over time, since the islets are exposed to the combined stress of high glucose and IL-1β. This is not explicitly incorporated in our model; S2 is set to a certain constant value at the beginning of our simulations but could be made a function of exposure time to fit the experimental transitions.The effect of IL-1 receptor antagonist on islet fateBecause IL-1β is an established part of T2D pathogenesis, different types of treatments targeting the cytokine have previously been investigated. Clinical studies with anakinra, a recombinant IL-1 receptor antagonist (IL-1Ra) have showed improvements in glycemia and β-cell function22. Remarkably, a 3-month anakinra treatment had a durable effect, lasting for over 39 weeks after anakinra withdrawal. Other studies targeted the increased level of cytokines through IL-1β antibodies, with a similar improvement of β-cell function and lowered inflammation39,40. While the mechanism behind the persistent effects after drug withdrawal is unknown, it has been hypothesized that it could result from interrupting the auto-inflammatory positive feedback loops41. We set out to test this hypothesis with our model. The result of blocking IL-receptors with the antagonist effectively corresponds to increasing the activation threshold of NF-κB in the model. This can be done in two ways: preventive treatment, where the drug is given before islets have transitioned into the locked state or reactive treatment, where the drug is given after the islets transitioned into the locked state.We find that preventive treatment offsets the transition by shifting the locked state to higher S2, Fig. 4A. In other words, islets can withstand higher levels of glycemia or free fatty acids in the presence of the antagonist. To assess the degree of dose-dependency of the simulated treatment, we have normalized the results in Fig. 4A by the untreated results (corresponding to 1.0 KI), see Fig. 4C. The normalized values can be thought of as islets “rescued” by the treatment. We find that the effect is dose-dependent in several ways: first, within a certain range of S2 (log2S2 = −0.92:−0.9), the fraction of “rescued” islets increases with increasing dose. Second, the range of S2 where the treatment has an effect also increases with the dose. While the trend is the same for the case of reactive treatment, where antagonist is added after islets have already locked, we find that reactive treatment is less efficient. Rescue of about 50% of the islets required 10-fold higher doses (higher activation threshold) compared with the preventive treatment.Fig. 4: Simulated effect of preventive and reactive treatment with IL-1β antagonists.A The transition of islets from transient to persistent inflammation with increasing levels of S2. Preventive treatment is simulated by changing the activation threshold, KI, prior to S1 exposure. Increasing KI offsets EC50 to higher S2. Simulations of the untreated case (1KI, red). KI is increased by ΔKI 0.1% (yellow), 0.25% (green) and 0.5% (blue). B Reactive treatment is simulated by changing activation threshold after S1 exposure. The results of the simulations in the untreated case with (1KI, red) and increased activation threshold by 1% (yellow), 5% (green) and 7.5% (blue). C, D The fraction of “rescued” islets represented by the relative change (treatment − no treatment) in fraction of locked islets. The reactive treatment has a weaker effect than preventive treatment, and the range of \({{\rm{S}}}_{2}\) where the treatment has the maximum effect is dose-dependent. To identify EC50, all data was fitted by sigmoidal curves (colored lines), see “Methods”. Results are from simulations in islets of radius 6.7 cells.While interfering with the positive feedback loop through both preventive and reactive treatments can rescue the islets from locking, we find that removing the treatment (resetting the KI back to pre-treated values) would not have long-lasting effects as long as S1 and S2 are unchanged. It is however known that the IL-1R antagonists improve glycemia and FFA levels in T2D patients41. This corresponds to lowering S2 prior to treatment withdrawal in our model, and, in this case, the effect of single islet unlocking with IL-1RA would be long-lasting.Islet shape, size and internal geometry play crucial roles in the fate of isletsHuman islets are diverse in size and geometrical arrangement of β-cells. Small islets tend to have a core-mantle-shaped structure with a homogenous core of β-cells, surrounded by α- and δ cells. Larger islets are more complex, as the core-mantle geometry is replaced by a heterogeneous mixture of α-, β- and δ cells and fenestrated capillaries24. Notably, it is large islets that are the first to be lost in T2D patients24. While the mechanism is currently unknown, there is strong evidence that inflammation plays an important role. First, Ehses et al.15 found that in murine models of T2D, macrophages infiltrate large islets first. Second, Youm et al.25 showed that ablation of NLRP3 inflammasome in chronically obese mice protected large islets from inflammation-induced death.These observations, together with our findings, suggest that the larger islets may have a higher propensity to lock and, as a consequence, would be compromised at lower levels of of S2. In order to test this, we first varied the size of the core-mantle-shaped islets and monitored the fraction of persistently responding islets with increasing S2 as described in Fig. 2D. We found that the transition from transient to persistent IL-1 response indeed occurs at lower values S2 when islets are larger (Fig. 5A).Fig. 5: Islet size and shape influence its propensity to lock.A Larger islets transition to the persistent mode at lower S2. The fractions of persistent islets (simulation resutls) are color-coded according to the islet sizes, red, orange and yellow corresponding to islets consisting of 199, 163 and 127 cells. B Islets with few or no ISCs in the core sustain high values of S2 without locking. Here green and blue mark the simulated fractions of persistent islets for respectively “human” and a donut-shaped islets. The islet to the left is a cross-section of a human islet configuration experimentally reported in ref. 18. The results are shown for S1 = 25 h−1 at each source, with number of S1 sources ns = 10. C Large core-mantled shaped islets are also more prone to locking when exposed to increasing S1. S1 at each source is 25 h−1, log2(S2) = 0.5. The color correspond to the islet sizes and configurations as defined in A and B. D The sources of S1 positioned deep in the bulk of the ISCs increase islet propensity to lock. The frequency of S1 at a given position is color-coded in islets that either respond transiently (left) or persistently (right). Results are sampled over 200 simulations where the position of 10 sources of S1 was chosen randomly among ISCs while other parameters were kept constant (S1 = 25 h−1). E The density of ISCs in the vicinity of S1 sources is estimated by the average number of ISCs within diffusion length \(\lambda =\sqrt{D{\tau }_{I}}\) from S1. For all shapes, the number of ISCs next to S1 sources is significantly smaller in transiently responding islets. (Kolmogorov-Smirnov test, comparing the distributions for each islet fate, see “Methods”). All simulations were run with p = 1400.To further quantify how the transition depends on the shape of the islet, we have implemented an experimentally reported configuration from ref. 18 (Fig. 5B, left). Despite a relatively large number of cells (159 cells), this configuration is more robust and can sustain higher levels of S2 (log2(S2lock) = 0.75) without locking compared to smaller islets with core-mantle configuration (127 cells, log2(S2lock) = 0). Consistent with this result, islets with ISCs distributed in a donut-shape, transit to the locked state at even higher S2. The results are qualitatively the same, if we fix levels of S2 and increase S1, Fig. 5C. If we assume both signals are constant throughout the pancreas, larger islets would be at higher risk of entering a self-sustained inflammatory response, simply because the critical S2 level is lower, consistent with observations in T2D patients and murine models24. In our simulations, a greater mixture of endocrine cell types or capillaries that can act as sinks can contain this threat.For each islet shape, there is a parameter range where transiently and persistently responding islets coexist. This is not the result of S1 and S2 but rather the internal geometry of sources in the islet. The position of the external sources is random and is the only source of stochasticity in our simulations. Interestingly, the region of co-existence is much broader for the human- and donut-shaped islet, indicating that the position of external sources within these islet types has a greater impact on islet fate. Voids, corresponding to non-ISCs or capillaries, could be a factor that increases fate diversity since a central sink can create a delay in the excitation between ISCs at different ends of the islet. In compact islets, the excitation will spread more or less uniformly, because we use a high diffusion constant compared with islet size.To investigate why there is a co-existence of modes, we have quantified if some positions of S1 are more likely to induce persistent modes while others are transient. In Fig. 5D, we show the frequency of source positions in transient and persistent modes. We find that the local density of ISCs around the sources is higher in persistent compared to transient modes. This is quantified in the box plot of Fig. 5E, where there are significant differences in the average number of ISC neighbors within the diffusion length, \(\lambda =\sqrt{D{\tau }_{I}}\), between islets of different fates. Hence an inflammatory response located at the center of an islet is more critical than peripheral cell defects or a response close to capillaries, where the cytokines can be cleared from the islet.Robustness of the results to parameter choices and cell-to-cell heterogeneityIn order to test the robustness of the model, we conducted a parameter sensitivity analysis.As described above, the parameters kIN, kRN, kNR, and τR, have been hand-fitted to simulate a simplified model for transient activation of NF-κB with a peak around 2 h34. We focused on less constrained parameters related to the production, degradation and diffusion of IL-1β (Fig. S6). As the transition from transient to persistent inflammation in islets is one of the key results, we have used it as a main readout for assessing how sensitive the model outcomes are to parameter changes. For a given parameter set and given number of S1 sources (ns), we simulate islet response to transient exposure to S1 signal. We use levels of IL-1β after S1 signal has been removed as a marker for the islet being in a locked (persistently elevated) or unlocked (returning to pre-stimulated level) state.To systematically quantify the transition point, we have fitted the outcomes of such an S1-scan to a step-like sigmodal function \(f\left(x\right)=V_{\max} \frac{{\left(x/K\right)}^{10}}{1\,+\,{\left(x/K\right)}^{10}}\). The fitted \(V_{\max}\) and \(K\), represent the levels of the IL-1 in locked islets and the transition point (number of S1 sources at which islet transitions to a locked state), respectively.We found that for the tested perturbations, the model still recapitulates the characteristic transition from transient (unlocked) to persistently elevated (locked) IL-1β levels in single islets and that the probability of “locking” increases with increasing level of Signal 1. In regards to sensitivity, parameters fall into two categories, the less sensitive τI, τI’, and D and one sensitive parameter, p. The less sensitive are the ones related to degradation and diffusion of IL-1 (different values of the parameters τI, τI’, and D). The variations in the transition point, K, do not seem to be sensitive to the exact values of these parameters (Fig. S6) but appear to be primarily dominated by the exact positioning of S1 sources. On the other hand, increasing p-values (describing the NF-κB mediated production of de novo IL-1 and effectively corresponding to the strength of the IL-1 positive feedback) correlated with increasing probability of locking (Fig. S6D). This is in line with our earlier findings27,29, where the strength of the positive feedback was identified as a bifurcation parameter. Extreme values of p lead to a disappearance of the co-existence regime: when p is small enough, no islets will lock, and when p is high enough, all islets will lock for the entire range of Signal 1 dose that was investigated. The level of IL-1 in the sustained state, Vmax, also increases with increasing p-values (Fig. S6E).The diffusion constant used in the main results was based on the estimate of diffusion of IL-1β from its molecular weight in cell-free medium; however, binding of IL-1 in the extracellular matrix is likely to result in a lower effective diffusion42. We accounted for this by checking that our results do not change qualitatively if the diffusion constant is further decreased up to 100-fold, see Supplementary Movies 7–9. Lowering the diffusion constant shifts and decreases the relevant parameter regimes (S1, S2 and p), since less cytokine is removed from the islets. This scenario can lead to interesting phenomena, as seen in excitable media, such as traveling waves or spiral patterns. In rare instances, source position can generate waves that travel around central arteries (see Supplementary Movies 10 and 11).Islets of Langerhans are subject to cell heterogeneity—not only in different cell types but also variations in protein levels inside and outside the cells. Noise in parameters such as inflammasome activation or clearance of inhibitor can have a large impact on the local dynamics, and we have therefore tested that the main results of our model are also valid when we introduce cell-to-cell heterogeneity in model parameters. In practice, we introduce such cell heterogeneity by modeling different parametrizations in each cell by drawing individual values of the parameters from a normal distribution with a mean given by the reference parameters set (Table 1 in “Methods”) and a width given by the coefficient of variation, CV. We choose random values for the parameters τI, τI’, p, and the activation threshold KI (reference value is equal to one, as we model the rescaled equations, Fig. 1C, also see “Methods” section). We do not choose random values for the diffusion constant, since this is inherently a global variable. Notice also that we have not performed a parameter scan of the activation threshold, since the model is scaled according to this parameter. It is however relevant to investigate cell-to-cell variation in the activation threshold, as studies have shown that cells do exhibit variationsThe results of simulations including cell-to-cell variability are shown in Supplementary Fig. S7. We found that the model still predicts transient and persistent production of IL-1β and that the co-existence regime widens when cell-to-cell heterogeneity is introduced (Supplementary Fig. S7C–E).
Tipping-point transition from transient to persistent inflammation in pancreatic islets
