Model simulations and parameter sampling were carried out with COPASI version 4.39 (Hoops et al., 2006; Bergmann et al., 2017, RRID:SCR_014260), Virtual Cell version 7.5.0 (Schaff et al., 1997; Moraru et al., 2008, RRID:SCR_007421), Tellurium version 2.2.7 (Choi et al., 2018) that uses libRoadRunner version 2.3.2 (Welsh et al., 2023, RRID:SCR_014763), and AMICI version 0.11.25 (Fröhlich et al., 2021), which was accessed through runBioSimulations (Shaikh et al., 2021, RRID:SCR_019110). The model file was constructed with python scripts using the BasiCO package that interfaces with COPASI (Bergmann, 2023). Simulations were run at the local high-performance computing cluster using the Cloud-COPASI web interface (Kent et al., 2012). Results were visualized with COPASI, with Gnuplot version 5.4.3 (Williams and Kelley, 2022, RRID:SCR_008619), or with the Python libraries Seaborn (Waskom, 2021) (RRID:SCR_018132) and Matplotlib (Hunter, 2007, RRID:SCR_008624). The SBGN diagram of Figure 1 was created using Cell Designer version 4.4 (Funahashi et al., 2003, RRID:SCR_007263) and then edited with Inkscape version 1.1 (RRID:SCR_014479).

The model used here is the segment polarity network model described by von Dassow et al. (2000). Briefly it represents a hexagonal array of cells, where each cell can express various genes (wingless, engrailed, hedgehog, cubitus interruptus, and patched) and where their protein products interact within a cell, and across neighboring cells. Figure 1 depicts the interaction network using the SBGN standard (Le Novère et al., 2009; Touré et al., 2021). Note that von Dassow et al. (2000) analyze two versions of this model, one having less interactions than the other. Here we only look at their full model (i.e., including the dashed arrows in the diagram of their Box 1). Since a 1 × 4 grid of cells is enough to replicate the results (von Dassow et al., 2000), that was used here to obtain all results.

My implementation of the model was first created for the widely used software COPASI (Hoops et al., 2006; Bergmann et al., 2017) through a Python script that creates a model with arbitrary number of cells at the user’s desire. A second script was created to generate the same model with only one cell, where the interacting species from neighboring cells are included as fixed concentrations. COPASI generates the full set of ODEs automatically based on the network and reaction kinetic rate laws. Unlike the SBML version from Sethna (2008), here we have the full reaction network, not just the differential equations. A small formal difference between this version and the original SPN model, is that COPASI expresses ODEs in terms of the species amounts rather than concentrations, but since the cell volumes are not variable this makes no difference and both sets of equations are equivalent.

The model makes extensive use of Hill-type functions where various terms appear in the form base ^exponent . This is often problematic in IEEE floating point since, for non-integer exponents, those operations are carried out based on the equivalence: base exponent = e exponent × log base . (1)

Therefore, calculations fail when base is negative, even if infinitesimally small (generates a NaN, which in COPASI is translated to an error “Invalid state”). Unfortunately, due to the nature of predictor-corrector integration algorithms, this can easily happen during a time course integration if one species concentration becomes very close to zero. In order to avoid this problem one can use a kind of “guarded” exponentiation: base exponent ≃ max ϵ , base exponent , ϵ > 0 . (2)

Applying this protection to the model changes the rate laws. For example, the rate law for transcription with inducer-repressor pair changes from the original: V ⋅ I ⋅ 1 − R h 2 k 2 h 2 + R h 2 h 1 k 1 h 1 + I ⋅ 1 − R h 2 k 2 h 2 + R h 2 h 1 (3) to the alternative: V ⋅ I ⋅ max ϵ , 1 − max ϵ , R h 2 k 2 h 2 + max ϵ , R h 2 h 1 k 1 h 1 + I ⋅ max ϵ , 1 − max ϵ , R h 2 k 2 h 2 + max ϵ , R h 2 h 1 . (4)

The terms k 1 h 1 and k 2 h 2 are not protected by a “guard” because k ₁ and k ₂ are constants that are always positive. In the results presented here I have used ϵ = 10^–80, which reduced the incidence of simulations with NaNs from around 10%–0.1%. von Dassow et al. (2000) did not describe how they avoided this problem within the software Ingeneue. Use of these alternative rate laws was necessary for the random parameter sampling, but for specific time course simulations one can almost always use the original rate laws as described in von Dassow et al. (2000).

Several aspects of the original SPN model were not fully described by von Dassow et al. (2000) and I have had to resort to later publications to infer what they could be. For the sake of complete transparency, here are all the details that had to be inferred from sources other than the original article.

• Parameter H _EWG does not feature in the differential equations of the Supplementary Material S1 or in von Dassow and Odell (2002), instead there the proteins EWG and IWG have the same half-life (H _IWG ). However the parameter is clearly described as one of the 48 parameters sampled in Meir et al. (2002), from the same group. Thus in my implementation EWG has its own half-life H _EWG .

As described in the Interoperability section of Results, below, the model can be exported from COPASI in standard formats, particularly the systems biology markup language (SBML, Hucka et al., 2003; Keating et al., 2020) and the OMEX format (Bergmann et al., 2014) containing a SBML file for the model and a SED-ML (Waltemath et al., 2011b) file with the simulation specification.

It is rather unfortunate that the term “reproducibility” has itself been used with various different meanings. This confusion in terminology was discussed in detail by Goodman et al. (2016), Plesser (2018), Miłkowski et al. (2018), and especially Barba (2018). As previously (Mendes, 2018), I will follow the definitions of Goodman et al. (2016), which specifies three distinct types of reproducibility.

• reproducibility of methods requires one to be able to exactly reproduce the results using the same methods on the same data;

Because the software Ingeneue, originally used to build and simulate the SPN model, has now disappeared from circulation, reproducibility of methods can no longer be effectively carried out. In a later publication von Dassow and Odell (2002) appear to have reproduced the results with the same software (see Table 1), however since these are the original authors, that can hardly be seen as independent verification. Of all the works listed in Table 1, only Ingolia (2004) and Ma et al. (2006) can be seen as independent reproductions of the original results. Unfortunately those two publications used their own C programs but did not publish them. It was work in Sethna’s lab (Gutenkunst et al., 2007; Daniels et al., 2008) that resulted in an electronic version of the model being created in the SBML format that is still available (see notes to Table 1), and which was re-used by Marazzi et al. (2022). However this SBML implementation coded the ODEs directly without representing the reaction network, an important limitation.

I attempted to reproduce the results of Table 1 in von Dassow et al. (2000), displayed in our Table 2. Overall these results match the original ones fairly well. There are some discrepancies in two samplings, but these are likely due to the uncertainty on the actual initial values, as pointed out in Methods. Bear in mind that these are very small samples of a 48-dimensional parameter space and the differences may just be due to random sampling. Figure 2 displays the succesful parameter sets in the sampling with crisp initial conditions, corresponding to Figure 2A in von Dassow et al. (2000). Careful comparison between the Figure and the original one reveals similar distributions. For example, in both cases κ _CNen rarely takes large values. The conclusions taken by von Dassow et al. (2000) would not change if their Figure 2A was substituted by this Figure 2. Taking these results together, I propose that the current implementation of the SPN model matches the results of the original—reproducibility of results.

To demonstrate that this implementation of the SPN model is interoperable across different software, a specific time course was chosen to be run by several simulators (herafter named timecourse1). One of the successful parameter sets generated in the random sampling with the “degraded” initial condition was chosen and saved as a native COPASI file, an SBML Level 3 Version 1 file (Hucka et al., 2018), and an OMEX file (Bergmann et al., 2014). Both the COPASI and OMEX files include the specification of the time course (end time of 1100 time units, sampled every 5 time units), though the SBML file requires that time course to be specified separately in the destination simulator.

Timecourse1 was simulated in four different software tools: COPASI, Virtual Cell (Schaff et al., 1997; Moraru et al., 2008), Tellurium (Choi et al., 2018), and AMICI (Fröhlich et al., 2021). It was run locally with COPASI, Virtual Cell, and Tellurium, and through the web service runBioSimulations (Shaikh et al., 2021) with AMICI. COPASI used the native file format, Tellurium used the SBML (through a small Python script runTellurium.py), while Virtual Cell and AMICI used the OMEX file.

Figures 3, 4 display the time course simulations obtained with four different software. There are no visible differences in the trajectories displayed confirming that these packages are all equally able to reproduce the results. Note that different ODE solvers were used by each one: COPASI used LSODA (Petzold, 1983), Virtual Cell used a fixed-step size Adams-Moulton method (Han and Han, 2002), Tellurium used CVODE (using the Adams-Moulton variable order, variable step size method) and AMICI used CVODES, both part of the SUNDIALS suite (Hindmarsh et al., 2005).

To promote findability and accessibility, the model files and associated scripts are made available through the following channels: a) a GitHub repository https://github.com/pmendes/models/tree/main/vonDassow2000, b) a Zenodo accession DOI (doi:10.5281/zenodo.7772570), c) a submission to the Biomodels database (MODEL2304060001), and d) model files deposited in the database of public Virtual Cell models. Note that the complete result files are only accessible through Zenodo since several files were larger than the limit at GitHub.

To demonstrate how the model can be reused for different purposes, I decided to ask the question “how often do parameter sets of the SPN model have multiple steady states?” Earlier von Dassow and Odell (2002) and especially Ingolia (2004) proposed that the robustness of pattern formation in the SPN model is due to multi-stability of steady states. Ingolia (2004) showed this in SPN models of a single cell (where the interacting species from the neighboring cells are kept constant). Here I investigate the answer to this question in a 1 × 4 array of cells. The strategy I used is as follows.

1. Generate p random sets of parameter values;

COPASI can easily to carry out such a study directly with the Parameter scan and Steady state tasks. The steady state task was applied here disabling the Newton method and therefore only using ODE integration to find the steady state reachable from the initial conditions (the steady state resolution was set to 10^–4 and the criterion used was “distance and time”). With the parameter scan task, 5,000 random parameter sets were sampled, using the same rules as in Section 3.1 above. Then, for each parameter set, it sampled 15 random initial conditions. Since we use a model of 1 × 4 array of cells, the initial conditions are composed of 132 species concentrations that were sampled in the interval [0,1].

From the 5,000 random parameter sets generated, 3,387 had at least one steady state (the remainder are likely to contain limit cycles, but this was not investigated). Of those 3,387 parameter sets with steady states, 498 contained more than one steady state. This rate of 1/10 parameter sets displaying multistability is not entirely surprising given the study by Ingolia (2004) which highlighted the positive feedbacks contained in the SPN model. Nevertheless it is interesting to investigate if these 498 parameter sets have special characteristics versus the other 2,889 that have only one steady state.

The distributions of parameter values that support multiple steady states was compared with those that appear to only support a single steady state. Calculation of the relative change in the median values for each parameter in the single steady state set versus the multiple steady state set revealed that only κ _CNptc shows a large difference, with a median 5-fold larger in the multiple steady state set than in the single steady state set. Three others have much lower differences: κ _CNen 0.7-fold smaller, κ _CIptc 0.46-fold smaller, and HH ₀ 0.45-fold smaller. The other 44 parameters have smaller differences. Figure 5 depicts the distributions of values of κ _CNptc and κ _CNen for the two data sets. Supplementary Figures S1–S3 depict histograms for all of the 48 parameters. There seems to be very few parameter sets that lead to multiple steady states with low values of κ _CNptc , while many more have high values for this parameter. This suggests that in order to achieve multiple stability the repression of patched (ptc) transcription by the truncated protein product of cubitus interruptus (CN) should be weak. Note that there is another negative feedback loop between these two genes, through induction of ptc transcription by the full length cubitus interruptus protein (CI).