Exploring the mechanistic and temporal regulation of LRP6 endocytosis in canonical WNT signaling

WNT signaling regulates central developmental processes of the cell, including cell fate, cell proliferation, cell migration and adult homeostasis. At the same time, aberrant or deregulated forms of WNT signaling are involved in a number of human cancers and developmental disorders (Moon et al., 2004; Clevers and Nusse, 2012; Logan and Nusse, 2004). Among the different WNT signaling pathways, the canonical WNT or WNT/β-catenin signaling pathway is the most intensively studied, in vitro (MacDonald et al., 2009) as well as in silico (Lloyd-Lewis et al., 2013). Central to initiation of the pathway is the successful phosphorylation and accumulation of the pathway-specific receptor low-density lipoprotein receptor-related protein 6 (LRP6) and its subsequent endocytosis. The endocytosis of LRP6 is an important regulatory mechanism that allows the cell to decide between signal activation and attenuation in general, as well as to differentiate between canonical and non-canonical signaling cascades (Hagemann et al., 2014; Feng and Gao, 2015). Various studies have demonstrated that LRP6 endocytosis is required for assembly of LRP6 signalosome and LRP6 phosphorylation and therefore is vital for the accumulation and transcriptional activity of β-catenin (Liu et al., 2014; Demir et al., 2013; Blitzer and Nusse, 2006).

According to the current understanding of this process in vertebrates, the WNT/LRP6 complex is internalized in microdomains (lipid rafts) in a caveolin-dependent fashion (clathrin-independent endocytosis; CIE) after the binding of WNT (Yamamoto et al., 2006; Özhan et al., 2013; Demir et al., 2013). Once internalized, the caveolin-containing endosome couples with the β-catenin destruction complex in which LRP6 is phosphorylated (Jiang et al., 2012). Depending on the ligand or specific adapter proteins, LRP6 can also be subject to clathrin-mediated endocytosis (CME) (Yamamoto et al., 2008; Kagermeier-Schenk et al., 2011; Jiang et al., 2012). However, in this case, LRP6 is not phosphorylated and no interaction with the destruction complex occurs. Instead, LRP6 and WNT are primarily internalized for ligand/receptor clearance and signal attenuation. To summarize, according to the aforementioned studies, the activity state of canonical WNT signaling is primarily determined by the endocytic pathway employed after WNT–LRP6 binding. Clathrin-independent (i.e. caveolin-mediated) internalization of LRP6 promotes β-catenin activity (i.e. its accumulation and transcriptional activity in the nucleus) and when LRP6 is shifted to clathrin-dependent internalization, β-catenin activity is reduced or shut down completely (Kagermeier-Schenk et al., 2011; Sakane et al., 2010).

In contrast to this, another set of studies provides an opposite picture of endocytic mechanisms regulating WNT/β-catenin signaling in which CME instead of CIE is the dominant internalization pathway for the WNT/LRP6 complex. As described by Blitzer and Nusse (2006) and recently reviewed by Brunt and Scholpp (2018), WNT/LRP6 is primarily internalized through dynamin- and clathrin-dependent endocytosis (CME). Perturbing this endocytic process resulted in attenuated WNT/β-catenin signaling, and therefore provided evidence that it is necessary for substantial β-catenin accumulation. Thus, central questions regarding the exact mechanism of endocytosis, such as which internalization route is activated under which circumstances and conditions, are still under debate (Brunt and Scholpp, 2018; Feng and Gao, 2015; Gagliardi et al., 2008).

Here, we use computational modeling to evaluate and discuss current hypotheses regarding the internalization of LRP6. We put particular emphasis on the role of membrane microdomains (also termed lipid rafts) because they control receptor activation, aggregation and endocytosis in WNT signaling as well as in many other signaling and endocytosis pathways (Sorkin and von Zastrow, 2009; Diaz-Rohrer et al., 2014). The membrane is structured into liquid-ordered (lo) and liquid-disordered (ld) domains.

The lo domains are often considered as membrane microdomains or lipid rafts. CME mostly occurs in ld domains (i.e. outside of microdomains) whereas CIE processes are typically considered as lipid raft-mediated, also because caveolin is primarily associated to lo domains (Lajoie and Nabi, 2010; Sakane et al., 2010; Le Roy and Wrana, 2005). Note that we use the terms ‘microdomains’ and ‘lipid rafts’ synonymously for membrane domains that are in the liquid-ordered state.

We incorporate existing hypotheses into models of endocytosis of different complexity and compare the temporal LRP6 internalization rates obtained from our simulations with experimental data gathered from various resources and studies.

Once validated, the model prospectively allows assessment of the spatiotemporal regulation of receptor endocytosis under various conditions depending on species, tissue and cell type; this has important implications for subsequent recycling and/or degradation dynamics and, hence, for overall cell signaling dynamics.

Here, several model implementations describing LRP6 endocytosis on various levels of detail are described and evaluated. Each model is a successive extension of the previous one. The temporal internalization dynamics of the WNT/LRP6 complex predicted by the model is compared with in vitro measurements obtained by Yamamoto et al. (2006), Khan et al. (2007) and Liu et al. (2014). We considered a parameter configuration of the model as valid if the fraction of internalized LRP6 receptors corresponded to experimental values at all time points. If no valid parameter configuration for a certain mechanistic model was found, another detail, such as an additional mechanism, cellular process or assumption, was added to the model, which was then reevaluated. Thereby, we derive insights into which cellular mechanisms are inevitable for the endocytosis of LRP6, allowing better characterization of this crucial process in canonical WNT signaling.

First, we evaluate the qualitative characteristics of the models that were built in the course of the study. These models describe basic mechanisms such as ligand–receptor binding, various internalization pathways and signalosome formation. The models are introduced in a pathway-independent manner to be applicable to different signaling pathways. Therefore, we start with a basic internalization model, which is successively extended. After this, these abstract models are parameterized for the WNT signaling pathway and the corresponding simulation experiments are presented.

For many ligand/receptor systems, endocytosis models can be implemented in a rather straightforward way and described by only a few reactions. A prominent example is endocytosis of the epidermal growth factor receptor (EGFR), whose internalization dynamics have been successfully described by a system of only two reactions: the formation of a ligand/receptor complex by reversible ligand–receptor binding and subsequent internalization of the receptor/ligand complex (Lauffenburger and Linderman, 1996). The dynamics of these two reactions are controlled by the association and dissociation rate constants k_a and k_d and the internalization rate constant k_e (Fig. 1A).

Note that in this model (A1) it is only possible to employ one endocytosis pathway: either CME or CIE. Because both endocytosis pathways differ significantly in their internalization dynamics, the value of k_e defines which endocytosis pathway is used in the given system (Goh and Sorkin, 2013). CME occurs at 0.2–0.4 min⁻¹, whereas CIE is typically in the range of 0.05–0.1 min⁻¹. Accordingly, a highly reduced model that solely describes the binding of a ligand (L) and a receptor (R), as well as the subsequent internalization of the resulting ligand/receptor complex (LR) can serve as a general (pathway-independent) description of ligand-mediated receptor internalization (Fig. 1A). Despite its simplicity, such a model, when parameterized accordingly, reproduces experimental data from pathways that employ either clathrin or caveolin/raft-mediated endocytosis as shown, for example, for EGFR internalization (Lauffenburger and Linderman, 1996).

Next, we consider a compartmental model that, as before, takes receptors and ligands, the cell and the extracellular environment into account. The model also structures the cell into cytoplasm and membrane. Thus, a compartmental hierarchy results (see Fig. 1B). Additionally, the membrane further separates into lipid raft and non-raft domains. This refers to the structural organization of the membrane into liquid-disordered (ld) and liquid-ordered (lo) domains. The lo domains are often considered as membrane microdomains or lipid rafts and are characterized by a local aggregation of saturated phospho- and sphingolipids and cholesterol (Lingwood and Simons, 2010; Sezgin et al., 2017).

In our model, membrane-integral receptors (R) can shuttle between raft and non-raft domains by means of lateral diffusion. To calculate the shuttling rate, we consider the membrane as a two-dimensional layer with lo domains being immobile circular-shaped entities within the membrane whose radius and coverage control the rate of receptor–domain association. We consider R to be (initially) homogeneously distributed across the membrane and to have no specific affinity to microdomains. The corresponding shuttling rate k_sh is solely based on diffusional association and can be approximated as described by Howard et al. (2001). As before, receptors can only be internalized when they form a complex with a ligand.

Typically, the internalization pathway depends on the localization of the ligand/receptor complex. This means that when located in non-raft domains, the ligand/receptor complexes are internalized via a clathrin-dependent pathway, whereas in raft domains, endocytosis is induced via the clathrin-independent caveolin pathway.

To account for this in model A2, we introduce two separate internalization reactions with corresponding reaction rate constants k_{e_nonraft} and k_{e_raft}. Thereby, it is also possible that both CME and CIE pathways are active at the same time as, for instance, observed by Di Guglielmo et al. (2003) (Fig. 1B).

Here, the compartmental model (as depicted in Fig. 1B) is coupled with a simple intracellular signaling model that includes the recruitment of cytosolic proteins, receptor activation and formation of signalosomes. Signalosomes are multifunctional protein complexes that contain membrane-localized receptors and cytosolic protein complexes. To include the most basic representation of signalosome formation in our abstract model, two additional reactions are added: the kinase (K)-mediated activation (phosphorylation) of the ligand/receptor complex (LR→LR_p) and the subsequent binding of cytosolic adaptor proteins (A) to the activated receptor complex (LR_p→C), as displayed in Fig. 1C. The establishment of signalosomes attenuates the diffusion and shuttling out of the microdomains. The assumption that the signalosome is internalized by the caveolin/raft-dependent pathway (by CIE within microdomains) also applies for A3.

Assuming that WNT-mediated endocytosis of LRP6 during canonical WNT signaling is regulated in a similar manner as endocytosis of EGFR, we transferred the previously described abstract models to the LRP6 internalization process during canonical WNT signaling. Therefore, the key players and mechanisms of the WNT-mediated endocytosis of LRP6 were mapped to the variables and reactions of the abstract models (see Fig. 2). The corresponding parameter values are listed in Table 1. The source code of the model implementations is accessible in the Github repository (see Materials and Methods). In the following sections, we describe and discuss pathway-specific assumptions and the challenges we encountered during the translation, as well as model configurations and the results of simulation experiments.

To arrive at a model implementation of A1 in the context of canonical WNT signaling (i.e. WNT model M1), we incorporate the following assumptions: (i) the extracellular ligand (L) and membrane-bound receptor (R) correspond to WNT proteins and LRP6 receptors; (ii) the interaction between LRP6 and frizzled (FZ) is modeled implicitly (i.e. we assume that LRP6 and FZ receptors are already in close proximity upon WNT binding and immediately form the WNT/FZ/LRP6 trimeric complex); (iii) the corresponding ligand/receptor complex (LR) is a simplified abstract representation of the WNT/FZ/LRP6 complex that can be internalized without regard for its phosphorylation state or association of additional binding partners such as axin or dishevelled (Dvl).

Importantly, according to abstract model A1, we focus only on the interaction between LRP6 and WNT, and the subsequent internalization of the receptor complex in M1.

Despite its crucial role in controlling canonical and non-canonical WNT signaling, we do not explicitly consider the presence of FZ for LRP6 receptor internalization in the context of canonical WNT signaling. Although several studies indicate that LRP6 is internalized in the form of a ternary complex comprising WNT, FZ and LRP6, the binding affinity of FZ to WNT is similar to that for LRP6 and both receptors are homogeneously distributed throughout the membrane (Sezgin et al., 2017; Bourhis et al., 2010; Brennan et al., 2004; Semënov et al., 2001).

Thus, the impact of FZ on LRP6 internalization dynamics can be neglected, as both receptors are activated through WNT3a and internalized in a similar way and with similar dynamics (Yamamoto et al., 2006). Therefore, we assume that FZ is either implicitly part of the WNT/LRP6 complex or that WNT–LRP6 interaction alone is sufficient for inducing LRP6 internalization (Brennan et al., 2004). This means that we consider the binding of extracellular WNT to the membrane-integral LRP6 receptors and subsequent internalization of the resulting WNT/LRP6 complex in terms of two reactions. For these rate constants, we apply WNT pathway-specific parameters (i.e. the association and dissociation rate constants for WNT–LRP6 interaction; Bourhis et al., 2010) and use either CME- or CIE-specific values for the internalization rate constant (Goh and Sorkin, 2013).

With this model (M1), we performed two separate simulation experiments: one simulation experiment with k_e values sampled from the parameter range of CME (0.2–0.4 min⁻¹) and one with k_e values sampled from the range of CIE (0.05–0.1 min⁻¹). Thereby, we recorded the number of LRP6 receptors (in complex with WNT) that were internalized after the indicated amount of time. The simulation results are depicted in Fig. 3A,B, together with the values that were obtained in vitro.

Neither in the range of CME nor in the range of CIE, were we able to find a model/parameter configuration that fit the experimental data derived from vertebrate cells (Yamamoto et al., 2006; Khan et al., 2007; Liu et al., 2014). In fact, both endocytic pathways were characterized by faster temporal dynamics than observed in vitro. Because the parameters employed in the model (i.e. association and dissociation constants) corresponded to experimental measurements from previous studies, this finding indicates the involvement of further mechanisms that restrict the availability of LRP6 receptors for endocytosis.

WNT signaling and corresponding endocytosis depend on the distribution of LRP6 in lo and ld membrane domains (lipid raft and non-raft fractions) (Sezgin et al., 2017; Sakane et al., 2010; Kagermeier-Schenk et al., 2011). Even though LRP6 is not specifically enriched in lipid raft fractions, its localization has a crucial impact on the internalization dynamics (Yamamoto et al., 2006; Sakane et al., 2010; Özhan et al., 2013). Accordingly, the structural organization of the membrane into lo and ld domains and the corresponding receptor distribution is reflected in M2 and adopted for the WNT pathway-specific implementation M2. Although the assumptions made in M1 also apply for M2, the following assumptions are defined in addition to M1: (i) approximately 30% of the membrane is occupied by lo membrane domains (lipid rafts); (ii) because LRP6 is homogeneously distributed in the membrane, on average 30% of total LPR6 receptors are located in lipid rafts; (iii) as in the abstract model A2, the internalization route depends on the localization of the WNT/LRP6 receptor complex.

Accordingly, we performed simulations with either one or both internalization pathways active. In the latter case, the preferred internalization route of the receptor is solely dependent on its distribution to raft and non-raft domains. Still, we were not able to reproduce the experimental measurements (Fig. 4A–D). However, we observed a significant difference to the previous (basic) model M1 when we considered the model version in which only CIE in lo domains is active (Fig. 4B). Here, the predicted fraction of internalized LRP6 receptors was lower than the experimentally derived values. Because CIE is now restricted to receptors located in raft domains, only a fraction of LRP6 receptors (∼30%) is subject to caveolin-mediated internalization. This attenuates the internalization dynamics, in particular during the first minutes after stimulation. In contrast to that, both model versions in which CME is active, either alone (in non-raft fraction) or in combination with CIE, yielded much higher temporal internalization dynamics than observed in vitro (Fig. 4A; Fig. S1A).

Based on these results, we hypothesized that a model in which CIE can occur outside of lipid rafts would lead to a better fit of the experimental data (Fig. 4C,D). This assumption is based on several studies (Johannes et al., 2015; Hansen and Nichols, 2009) that reported the occurrence of CIE in non-raft domains, hence independent of lipid rafts. In contrast to this, clathrin-dependent internalization of raft-associated LRP6, or raft-associated receptor complexes in general, has not yet been observed. Therefore, we consider that clathrin-dependent internalization of raft fractions is highly unlikely and did not consider this model configuration in the simulation. Accordingly, we applied CIE-based parameter values for the internalization rate of WNT/LRP6 complexes in both membrane domains (i.e. inside as well as outside of lipid rafts; Fig. 4C,D). We performed simulations with the same three scenarios as considered before: active internalization in either ld or lo domains or concurrently in both domains. The simulation results met either the temporal dynamics of the early onset of internalization (10–20 min after stimulation) (Fig. 4C,D) or the later time points (60 min after stimulation) (Fig. 4B), but not both of them. Because no parameterization of the described compartment-based model delivered correct or realistic internalization dynamics of the WNT/LRP6 complex, we conclude that an additional delaying process or mechanism needs to be considered in the process of endocytosis.

A candidate for the aforementioned delaying process is the formation of the signalosome [i.e. the recruitment and binding of cytosolic proteins such as casein kinase 1 (CK1), Dvl or axin to the receptor complex] and its subsequent activation (Bilic et al., 2007), as described in A3. To account for these processes in our pathway-specific WNT model (M3), we combined the compartmental membrane model described above with our previously published and validated model of canonical WNT/β-catenin signaling (Haack et al., 2015). Thereby, central membrane-related processes of canonical WNT signaling, such as raft-dependent phosphorylation and accumulation of LRP6, the recruitment of axin and subsequent formation of signalosomes in raft domains, are integrated into intracellular β-catenin signaling processes. Additional assumptions and simplifications are as follows: (i) Axin, a representative member of the destruction complex, is recruited to the membrane and binds to the phosphorylated LRP6 receptor complex. Thus, the LRP6 signalosome is represented by a simplified form comprising phosphorylated LRP6, WNT and axin. (ii) Although the WNT/LRP6 complex can shuttle between raft and non-raft domains, the LRP6 signalosome cannot shuttle between membrane domains. (iii) Similar to the previous internalization model, the internalization route is determined by the domain association of LRP6.

Note that raft-dependent endocytosis requires successful signalosome formation (i.e. the LRP6 signalosome complex is internalized only after LRP6 has been phosphorylated and axin has bound to the phosphorylated LRP6 receptor complex). Again, we ran simulations with all possible combinations of internalization mechanisms (despite clathrin-dependent internalization in rafts). Model configurations with CME and CIE being constrained to the outside of rafts (Fig. 5A,B) yielded slightly faster and slower internalization dynamics, respectively. Additionally, neither process in which the clathrin-independent internalization was constrained solely outside (Fig. 5B) or inside (Fig. S1B) of rafts alone was able to reproduce the in vitro internalization dynamics.

Finally, we ran simulations in which CIE occurs within the raft domains (restricted to the signalosome) and where CIE can also occur outside of raft domains (Fig. 5C). Intriguingly, we now have a model configuration that perfectly fits the experimental data. The best fit was achieved for the upper boundary values of ld internalization rate constants (k_{e_nonraft}) (i.e. for parameter values ranging between 0.075 and 0.1 min⁻¹). The model is characterized by CIE of the WNT/LRP6 complex outside of lipid rafts and CIE of the established LRP6 signalosome inside lipid rafts.

To ensure that our model is not overfitted to the given experimental data, we evaluated whether the model is able to account for the mechanisms of other endocytosis processes related to the canonical WNT signaling pathway. For this, we use the temporal dynamics of Dkk1-induced LRP6 endocytosis as measured by Sakane et al. (2010). Dkk1 shifts the endocytosis of LRP6 from clathrin-independent to clathrin-based internalization and, therewith, also changes the temporal dynamics of this process. Accordingly, we replaced WNT by Dkk1. The resulting model includes Dkk1/LRP6-specific association and dissociation rates (instead of WNT/LRP6), as obtained by Bourhis et al. (2010). CME is set as the only internalization pathway available for Dkk1/LRP6 complexes. Thus, we sampled from the parameter range of CME (0.2–0.4 min⁻¹) (Fig. 5D). The model immediately reproduced (i.e. without any additional changes) the internalization dynamics of LRP6 after binding to Dkk1 (Fig. 5D), as observed by Sakane et al. (2010).

Here, we employ a mechanistic modeling approach to explore and evaluate the cellular mechanisms required for LRP6 internalization. Starting with a relatively simple model that is controlled by only a small number of essential dynamics, we successively extended the model by cellular mechanisms until it satisfied a certain property. In our case, the desired property relates to the time series (trajectory) of LRP6 surface receptor concentrations in response to a WNT stimulus that corresponds to quantitative in vitro measurements (Yamamoto et al., 2006; Khan et al., 2007; Liu et al., 2014). However, the main question for any modeling and simulation-based study is how valid and explanatory the developed model is and how reasonable and meaningful the results derived from the corresponding simulation experiments are. The benefit of our WNT/LRP6 internalization models is that the kinetic rate constants of the essential dynamics are all defined by previous in vitro studies. The k_a and k_d values for the binding complex of LRP6 and WNT were determined by Bourhis et al. (2010) (see Table 1) and internalization rate constants for CME and CIE were obtained in several studies, as listed by Goh and Sorkin (2013). The same applies to the model extensions. The number of receptors as well as their distribution between membrane domains were experimentally derived by Bafico et al. (2001) and by Yamamoto et al. (2006) and Sakane et al. (2010), respectively. All intracellular dynamics of canonical WNT signaling are either based on measurements or have been thoroughly (cross-)validated after fitting by Haack et al. (2015).

For the time-dependent internalization dynamics of the surface WNT/LRP6 complex, we used Yamamoto et al. (2006), Khan et al. (2007) and Liu et al. (2014) for reference experimental data obtained from various (vertebrate) cell lines. In all cases, biotinylation assays were used. Note that we require quantitative temporal dynamics to fit our models because qualitative data (e.g. achieved by knockdown of either endocytosis pathway) does not provide enough information to uncover or explain underlying dynamics incorporated into mechanistic models. The small number of undefined rate constants in our models reduces the degree of freedom and allows us to test model configurations in a small range of parameter values. If no valid parameter configuration in the given parameter boundaries can be found for a certain model, the model itself becomes highly unlikely. Thereby, we challenge theories about the endocytosis of LRP6 during canonical WNT signaling and, with it, underlying hypotheses about the involved regulatory mechanisms. Several model extensions and specific assumptions were required to fit LRP6 internalization dynamics to experimental data, implying that LRP6 endocytosis is highly regulated and depends on a number of intertwined regulatory processes. These essential processes are highlighted in terms of pathway-independent abstract models (A1–A3).

After parameterizing the abstract models with WNT/LRP6-specific reaction rate constants (M1-M3, Fig. 2), we performed simulations with both clathrin-mediated (CME) and clathrin-independent (CIE) endocytosis pathways and compared the simulation results (i.e. the fraction of receptors internalized after an indicated amount of time) with experimental time-resolved measurements from various sources (Yamamoto et al., 2006; Khan et al., 2007; Liu et al., 2014). However, experimental data derived from vertebrate cells could only be reproduced by combining a model of central intracellular processes of canonical WNT signaling with a compartment-based membrane model in which we have representations of raft and non-raft domains (Haack et al., 2015). Because CME and CIE can occur in either of these domains (i.e. in ld and lo domains, respectively), such a compartment-based model allows the simulation of internalization of the WNT/LRP6 receptor complexes through both endocytosis pathways concurrently. Furthermore, integration of intracellular components of the WNT signaling pathway accounts for the delay caused by recruitment of members of the destruction complex and subsequent signalosome formation in lipid rafts before internalization (model M3, see Fig. 2C). Note that the compartmental model without intracellular processes (i.e. model M2) fails to reproduce experimental results under any valid parameter configuration. This indicates that the signalosome formation occurring exclusively in raft domains needs to be completed before internalization and that the resulting delay is a crucial part of the internalization dynamics of LRP6. This is in line with recent studies (Sezgin et al., 2017; Demir et al., 2013; Kim et al., 2013).

Intriguingly, we were able to reproduce the experimental results only under the assumption that internalization of WNT/LRP6 complexes occurs concurrently inside and outside of lipid rafts and is, in both cases, clathrin independent. At the same time, our model correctly predicts the DKK1-mediated shift from caveolin- to clathrin-dependent internalization of LRP6. This indicates that the internalization dynamics of our model are not overfitted to the internalization dynamics of WNT/LRP6 complexes, but allow the reproduction of other in vitro measurements of internalization processes related to LRP6 and WNT signaling.

Our findings thus confirm that, under normal conditions, LRP6 is subject to caveolin-mediated (clathrin-independent) internalization; however, according to our model, it is highly unlikely that caveolin-mediated LRP6 internalization occurs in microdomains or caveolae alone. Instead, most WNT/LRP6 receptor complexes are internalized directly after WNT binding through a CIE pathway outside of rafts. This means that our simulation results indicate a raft- and clathrin-independent endocytosis pathway with similar internalization dynamics as raft-dependent caveolin-mediated internalization. However, caveolin is raft associated, which makes it highly unlikely that caveolin also mediates the internalization of non-raft associated LRP6 (Azbazdar et al., 2019; Sonnino and Prinetti, 2009). Recently a flotillin-dependent internalization mechanism of LRP6 has been discussed (Yamamoto et al., 2017). Also, other potential clathrin-independent processes capable of building endocytotic pits have been discussed (Johannes et al., 2015), such as the clathrin-independent carrier (CLIC)–GPI-anchored protein-enriched early endosomal compartment (GEEC) pathway (Sabharanjak et al., 2002; Kirkham et al., 2005) or the formation of clathrin-independent pits induced by receptor clustering and mediated by short actin structures (Rao and Mayor, 2014). Clearly, this interesting aspect requires further investigation.

When comparing endocytosis between vertebrates and invertebrates, one has to keep in mind that invertebrates are not capable of producing or synthesizing sterols such as cholesterol and therefore have completely different membrane lipid compositions. In WNT signaling, this has crucial implications because cholesterol is one of the main components of membrane microdomains, which are crucial for activation and internalization of LRP6 in vertebrates. Without microdomains, it is quite likely that LRP6 is internalized by CME in invertebrates whereas, at the same time, LRP6 internalization in vertebrates might be mediated through lipid rafts without causing any contradiction.

Therefore, membrane lipid composition is another major aspect in signaling and, in particular, WNT signaling. Membrane composition determines how strongly the membrane separates into liquid-ordered and liquid-disordered phases as well as the properties of the resulting microdomains and how they affect receptor diffusion and, hence, receptor distribution.

Our model analyses further corroborate the importance of membrane composition in the context of WNT/β-catenin signaling and signal transduction in general.

This section briefly summarizes the methods used for developing the simulation models, for specifying and executing simulation experiments and for conducting the simulation study.

All simulation models (M1-M4) have been defined in ML-Rules, a multilevel rule-based modeling language (Maus et al., 2011). Rule-based modeling approaches have been introduced to allow more succinct and convenient modeling of biochemical systems (Danos and Laneve, 2004; Hlavacek et al., 2006; Faeder, 2011). Rule-based approaches have been and are frequently applied for simulating and analyzing signaling pathways (Salazar-Cavazos et al., 2020; Boutillier et al., 2018). The focus of ML-Rules on multilevel modeling and compartmental dynamics, such as the fusion and fission of compartments or endocytosis, has influenced language design. Species in ML-Rules can have attributes whose values can be of any type, such as integer or real. Attribute values can be accessed and updated by arbitrary (user-defined) functions. Similarly, the content of a compartment and the overall kinetics can be defined by these arbitrary functions.

The implemented simulation models make use of rule schemata and nesting. For example, in our coupled compartment-based models the central component of the signalosome, LRP6, is attributed with two different attributes: phosphorylation state and binding state. Individual LRP6 receptors diffuse between raft and non-raft domains irrespective of their phosphorylation or binding status.

Note that the models developed in this study are rather simple and we do not exhaust the entire functionality of ML-Rules. We chose ML-Rules as it allows us to easily extend the models (e.g. to model the growth of lipid rafts, and to include more detailed models of endosome dynamics, protein sorting and receptor recycling).

For a detailed presentation of ML-Rules, its syntax and semantics as well as a comparison with other rule-based approaches, the interested reader is referred to Helms et al. (2017).

The semantics of ML-Rules is based on continuous-time Markov chains (CTMC) (Helms et al., 2017). For ML-Rules, several simulation algorithms exist. For all simulation experiments performed in this study, we used a hybrid simulator (Helms et al., 2017). In these cases, only numerical integration methods were used. To confirm the applicability of this deterministic approach, we performed several analyses and test runs with all models to exclude the possibility of stochastic effects. For example, we checked whether all species are sufficiently abundant at any time and ascertained for each model that individual and mean trajectories of simulation runs executed with the Stochastic Simulation Algorithm (SSA) provide the same results as with the hybrid simulator. See Table 1 for an overview of the parameters used in the models. All simulation models are included in our GitHub repository, https://github.com/SFB-ELAINE/SI_LRP6_Endocytosis_Model.

All simulation experiments were executed in SESSL (Simulation Experiment Specification on a Scala Layer) (Ewald and Uhrmacher, 2014). SESSL provides a simulation-system-agnostic layer between user and simulation system, which can work with different simulation tools. It has been developed to support the compact specification and execution of a wide variety of simulation experiments. Proficient users can directly add missing features ‘on-the-fly’. As a domain-specific language embedded in the programming language Scala (https://www.scala-lang.org), SESSL offers extension points where user-supplied functionality can be included. In addition to simulation-based multi-objective optimization and statistical- and simulation-based model checking, different experimental design methods such as full factorial design, Latin Hypercube sampling or central composite design, and replication criteria such as maximal relative width of a confidence interval of simulation output, SESSL also supports sensitivity analysis (Warnke and Uhrmacher, 2018). Based on the bindings between SESSL and ML-Rules (Warnke et al., 2018), various simulation experiments were executed. During review, the editorial office asked to move any code into github repository. All data and scripts to replicate the plots are available in our GitHub repository, https://github.com/SFB-ELAINE/SI_LRP6_Endocytosis_Model. The files in the repository include the experimental data in data/; all experiment specifications in the subdirectories of experiments/; the R script ExcecuteSesslandPlotResults.R in the directory experiments/ to run all simulation experiments and create the plots that were used in the paper; and all models in the directory models/.

The documented results in Figs 3, 4 and 5 refer to the previously described parameter scans based on deterministic model simulations. To test whether the different models can be fitted to the in vitro measurements (Yamamoto et al., 2006; Khan et al., 2007; Liu et al., 2014) (Table S1), parameter scan experiments with parameter values for CME (0.2-0.4 min⁻¹) and CIE (0.05-0.1 min⁻¹) (Goh and Sorkin, 2013) were sampled for these model configurations. Note that in contrast to the endocytosis screens with purified WNT3a, the exact concentration of Dkk1 treatment underlying the experimental measurements obtained by Sakane et al. (2010) (Table S2) and displayed in Fig. 5D had to be fitted. This is because endocytosis screens with Dkk1 were performed with Dkk1-FLAG conditioned media instead of purified Dkk1.

We used the amounts of cell-surface LRP6 obtained after 30 min of treatment with different concentrations of purified Dkk1, as displayed in Fig. 3C of Sakane et al. (2010), to fit this parameter.

The different simulation models (M1-M4) used to explore cellular mechanisms of LRP6 internalization were developed by successive extension and revision. These processes are illustrated in Fig. 6, adopting the PROV-DM standard notation for provenance (Belhajjame et al., 2013). Capturing provenance (i.e. information about how a product has been generated) helps to assess the products of simulation studies, such as simulation data or simulation models (Ruscheinski and Uhrmacher, 2017). Fig. 6 gives a bird's eye view of the overall process and distinguishes between model building activities (solid rectangle) that generate simulation models (blue ellipse) and model analysis activities (dashed rectangle) that performed simulation experiments to check whether in vitro measurements can be reproduced. Therefore, the latter use in vitro measurements for comparing the simulation outputs (Yamamoto et al., 2006; Khan et al., 2007; Liu et al., 2014) and use in vitro measurements as inputs to scan the parameter values in a realistic range (Goh and Sorkin, 2013). It should be noted that the different roles of entities such as data or models can easily be included in provenance graphs (e.g. whether data have been used as inputs). This information was only withheld for clarity reasons. If provenance information is collected automatically, provenance graphs store detailed information about simulation studies as well as the involved activities and entities. These graphs can be queried to reveal interesting details on demand (Ruscheinski et al., 2019).

We would like to thank the reviewers for their valuable comments on the paper and Tom Warnke for his excellent technical support on simulation experiment specification in SESSL.

The peer review history is available online at https://jcs.biologists.org/lookup/doi/10.1242/jcs.243675.reviewer-comments.pdf