A mathematical model of p62-ubiquitin aggregates in autophagy.

Aggregation of ubiquitinated cargo by oligomers of the protein p62 is an important preparatory step in cellular autophagy. In this work a mathematical model for the dynamics of these heterogeneous aggregates in the form of a system of ordinary differential equations is derived and analyzed. Three different parameter regimes are identified, where either aggregates are unstable, or their size saturates at a finite value, or their size grows indefinitely as long as free particles are abundant. The boundaries of these regimes as well as the finite size in the second case can be computed explicitly. The growth in the third case (quadratic in time) can also be made explicit by formal asymptotic methods. In the absence of rigorous results the dynamic stability of these structures has been investigated by numerical simulations. A comparison with recent experimental results permits a partial parametrization of the model.

Introduction

Autophagy is an intracellular pathway, which targets damaged, surplus, and harmful cytoplasmic material for degradation. This is mediated by the sequestration of cytoplasmic cargo material within double membrane vesicles termed autophagosomes, which subsequently fuse with lysosomes wherein the cargo is hydrolyzed. Defects in autophagy result in various diseases including neurodegeneration, cancer, and uncontrolled infections (Levine and Kroemer 2008). The selectivity of autophagic processes is mediated by cargo receptors such as p62 (also known as SQSTM1), which link the cargo material to the nascent autophagosomal membrane (Danieli and Martens 2018). p62 is an oligomeric protein and mediates the selective degradation of ubiquitinated proteins. Its interaction with ubiquitin is mediated by its C-terminal UBA domain, while it attaches the cargo to the autophagosomal membrane due to its interaction with Atg8 family proteins such as LC3B, which decorate the membrane (Pankiv et al. 2007). Additionally, p62 serves to condensate ubiquitinated proteins into larger condensates or aggregates, which subsequently become targets for autophagy (Sun et al. 2018; Zaffagnini et al. 2018). It has been reported that this condensation reaction requires the ability of p62 to oligomerize and the presence of two or more ubiquitin chains on the substrates (Wurzer et al. 2015; Zaffagnini et al. 2018).

In this article a mathematical model for the condensation process is derived and analyzed. It is based on cross-linking of p62 oligomers by ubiquitinated substrate (Zaffagnini et al. 2018). A cross-linker is assumed to be able to connect two oligomers, where each oligomer has a number of binding sites corresponding to its size. As an approximation for the dynamics of large aggregates, a nonlinear system of ordinary differential equations is derived.

The oligomerization property of p62 has been shown to be necessary in the formation of aggregates (Zaffagnini et al. 2018): too small oligomers of Ubiquitin do not form aggregates (Wurzer et al. 2015).

The dynamics of protein aggregation has been studied by mathematical modelling for several decades, but most models consider the aggregation of only one type of protein, which gives rise to models belonging to the class of nucleation–coagulation–fragmentation equations, see e.g. Bishop and Ferrone (1984), Prigent et al. (2012) and Xue et al. (2008) for examples in the biophysical literature, and Collet et al. (2002), Laurençot and Mischler (2002), Banasiak and Lamb (2006) and Dubovskiĭ and Stewart (1996) for a sample of the mathematical literature. Contrary to these studies, the present work considers aggregates composed of two different types of particles with varying mixing ratios, which drastically increases the complexity of the problem.

In the following section the mathematical model is derived. It describes an aggregate by three numbers: the number of p62 oligomers, the number of cross-linkers bound to one oligomer, and the number of cross-linkers bound to two oligomers. The model considers an early stage of the aggregation process where the supply of free p62 oligomers and of free cross-linkers is not limiting. Since no other information about the composition of the aggregate is used, assumptions on the binding and unbinding rates are necessary. In the limit of large aggregates, whose details are presented in an appendix, the model takes the form of a system of three ordinary differential equations. Section 3 starts with a result on the well posedness of the model, and it is mainly devoted to a study of the long-time behaviour by a combination of analytical and numerical methods. Depending on the parameter values, three different regimes are identified, where either aggregates are unstable and completely dissolved, or their size tends to a limiting value, or they keep growing (as long as they do not run out of free oligomers and cross-linkers). In Sect. 4 we discuss the parametrization of the model and a comparison with data from Zaffagnini et al. (2018).

Presentation of the model

Discrete description of aggregates

We consider two types of basic particles: An aggregate is represented by a triplet , where i denotes the number of one-hand bound Ubi, j denotes the number of both-hand bound Ubi, and k denotes the number of p62. An aggregate will be assumed to contain at least two p62, i.e. , and enough both-hand bound Ubi to be connected, i.e. . Furthermore, an aggregate contains nk binding sites for Ubi, implying . A triplet satisfying the three inequalitieswill be called admissible. It is a rather drastic step to describe an aggregate only by these three numbers, since the same triplet might represent aggregates with various shapes. An example of an admissible triplet describing a unique aggregate topology is , representing a chain of p62. Adding one both-hand bound Ubi already creates a topological ambiguity: The triplet (0, k, k) can be realized by a circular aggregate or by an open chain, where one connection is doubled. Apart from missing topological information, the triplet (i, j, k) also lacks any information on the geometry of an aggregate. The reaction rates described below should be interpreted as mean values taken over all possible aggregate shapes described by a triplet (i, j, k). An extension of the model to include the shape information would require knowledge on the conformations of oligomers and cross-linkers, including chemical and mechanical properties, which is only partially available in the literature [see, e.g., Ciuffa et al. (2015) for p62 oligomers]. The modeling of aggregation reactions would then require molecular dynamics simulations, which would be rather complex even for small aggregates. For our goal of describing the dynamics of large aggregates such an approach seems to be prohibitively complex.

Oligomers of the protein p62, where we assume for simplicity that all oligomers contain the same number of molecules (see Sect. 4 for a discussion of this assumption). These oligomers are denoted by p62 and are assumed to possess n binding sites for ubiquitin each,

Cross-linkers in the form of ubiquitinated cargo, denoted by Ubi and assumed to have two ubiquitin ends each. When one end of a Ubi is bound to a p62, we call it one-hand bound, when both ends are bound we call it both-hand bound.

The reaction scheme

Basically there are only two types of reactions: binding and unbinding of Ubi to p62. However, depending on the situation these may have various effects on the aggregate, whence we distinguish between three binding and three unbinding reactions.

Addition of a free Ubi, requiring at least one free binding site, i.e. , (see Fig. 1): The reaction rate (number of reactions per time) is modeled by mass action kinetics for a second-order reaction with reaction constant and with the number [Ubi] of free Ubi. Since free Ubi and free p62 oligomers will be assumed abundant, their numbers [Ubi] and will be kept fixed and the abbreviation will be used (see Sect. 5 for a short discussion of the effect of lifting the abundancy assumption). This leads to a first-order reaction rate

Fig. 1

Examples for Reactions 1 (left) and 2 (right) with p62 in black, one-hand bound Ubi in green, two-hand bound Ubi in red, free particles in blue. Reaction 1: . Reaction 2: p62

Addition of a free p62, requiring at least one one-hand bound Ubi, i.e. : Analogously to above, we set and

Compactification of the aggregate by a Ubi binding its second hand, requiring at least one one-hand bound Ubi, i.e. , and at least one free binding site, i.e. : This is a second-order reaction with rate

Loss of a Ubi, requiring at least one one-handed Ubi, i.e. . This is the reverse reaction to 1: Its rate is modeled by

Loss of a p62 (leading to loss of the whole aggregate if ): This and the following reaction need some comments. They are actually both the same reaction, namely breaking of a cross-link, which we assume to occur with rate . However, this can have different consequences. Here we consider something close to the reverse of reaction 2. This means we assume that the broken cross-link has been the only connection of a p62 oligomer with the aggregate, such that the oligomer falls off. This requires the condition , meaning the possibility that the other binding sites of the lost oligomer are free of two-hand bound Ubi. It is not quite the reverse of reaction 2, since we have to consider the possibility that one-hand bound Ubi, , are bound to the lost oligomer. The conditional probability to be in this case, when a cross-link breaks, is zero for a very tightly connected aggregate where each oligomer is cross-linked at least twice, i.e. , and it is one for a very loose aggregate, i.e. a chain with . This motivates the model which is, as a function of , the linear interpolant between the values 1 and 0; notation: , leading to the rate In the framework of our model, should be a random number satisfying the restrictions where the upper bound should be obvious and the lower bound implies that the last condition in (1) is satisfied by the post-reaction state . We shall use the choice which can be interpreted as the rounded ( denotes the closest integer) expectation value for the number of one-hand bound Ubi on the lost oligomer in terms of the ratio between the number of available binding sites on the lost oligomer and the total number of available binding sites for one-hand bound Ubi in the whole aggregate. It is easily seen that in the relevant situation , i.e. , the choice (9) without the rounding satisfies the conditions (8). Since the bounds in (8) are integer, the same is true for the rounded version.

Note that we neglect the possibility to lose more than one oligomer by breaking a cross-link, i.e. the fragmentation of the aggregate into two smaller ones. On the one hand, this is a serious and actually questionable modeling assumption. Dropping it, on the other hand, would induce serious difficulties in the further development of our model. First, a decision would be required, which of the two post-fragmentation aggregates to follow. Second and more importantly, it would mean to allow for large jumps in the state of the aggregate, ruling out a differential equation model (as derived below), which is based on gradual state changes. An a posteriori justification of the no-fragmentation assumption will be provided by some of the results of the following section, showing that growing aggregates are tightly connected, making fragmentation very unlikely (Fig. 2).

Fig. 2

Examples for Reaction 3 (left, ), Reaction 5 (right, up, , ), and Reaction 6 (right, down, )

Loosening of the aggregate by breaking a cross-link, requiring at least one excess both-hand bound Ubi, i.e. : This is the reverse of reaction 3 with the rate which respects the requirement for a positive rate, because of

A deterministic model for large aggregates

The next step is the formulation of an evolution problem for a probability density on the set of admissible states (i, j, k). In this problem the discrete state is scaled by a typical value of and , assumed of the same order of magnitude. We also assume to be large, which means that we assume the concentrations [Ubi] and of free particles to be large enough (see Sect. 4 for a further discussion of these assumptions):This choice means that the orders of magnitude , , and of the forward reaction rates (2), (3), (4) are the same. It is then consistent with these definitions to introduce and to assume that takes moderate values. In the large aggregate limit , the new unknowns become continuous, and the equation for the probability density becomes a transport equation (see “Appendix A” for the details). It possesses deterministic solutions governed by the ODE initial value problemwhereis the limit of as . The conditions for admissible states are obtained in the limit of (1):implying, as expected,The equations satisfied by s and ,show that the conditions (14) are propagated by (12).

Analytic results

Global existence

Since the right hand sides of (12) contain quadratic nonlinearities, it seems possible that solutions blow up in finite time. On the other hand, the right hand sides are not well defined for . The essence of the following global existence result is that neither of these difficulties occurs.

Theorem 1

Let and . Let satisfy (14). Then problem (12) has a unique global solution satisfying as well as (14) for any . Also the following estimates hold for :

Proof

Local existence and uniqueness is a consequence of the Picard–Lindelöf theorem. Global existence will follow from the bounds stated in the theorem. Positivity of the solution components, of , and of is an immediate consequence of the form of the equations (12), (16), (17). This also implieswhich shows (18) by the Gronwall lemma. With (14), the equation for q in (12) impliesand another application of the Gronwall lemma and of (14) proves (19) and, thus, completes the proof of the theorem.

Long-time behaviour

The first step in the long-time analysis is the investigation of steady states. Although the right hand sides of (12) are not well defined for , the origin can be considered as a steady state sincehold for admissible states satisfying (14) and, thus, the right hand sides of (12) tend to zero for admissible . The following result shows that at most one other steady state is possible which, somewhat miraculously, can be computed explicitly.

Theorem 2

Let , , and letsatisfy . Then there exists an admissible steady state of (12) given bywith and . There exists no other steady state (besides the origin).

The origin is the only steady state with , since by (14), i.e. , implies . Assuming we introduceand rewrite the steady state equations in terms of and :From (24) we obtainwhich is substituted into the sum of (22) and (23):The option leads to , implying and, thus, , which contradicts (23). Therefore the second paranthesis has to vanish, leading to a quadratic equation for with the only positive solutionNow (24) implies the formula for stated in the theorem and we note that implies . We compute from by (25) and note that since . We then compute from the sum of (22) and (23):which proves . Finally we obtain the formula for from (23) as well as and .

For convenience below, the conditions in the theorem are made more explicit in terms of the parameters byThe steady state approaches the origin as , whereas all its components become unbounded as . This motivates the following.

Conjecture 1

With the notation of Theorem 2,

if , then all solutions of (12) converge to as ,

if , then all solutions of (12) converge to (0, 0, 0) as ,

if , then for all solutions of (12) we have as .

The conjecture has been supported by numerical simulations. Figures 3 and 4 show typical simulation results corresponding to the three cases. Partial rigorous results on the conjecture have been proven in the parallel work (Delacour et al. 2020). These include local stability of the zero steady state for . Since the right hand side of (12) lacks sufficient smoothness, this result cannot be proven by the standard linearization approach.

Fig. 3

Left: Convergence to the non-trivial steady state of Theorem 2. Simulation of an aggregate (p, q, r) of initial size (2, 4, 3) with parameters and , implying . Right: Instability of the aggregate. Simulation of an aggregate (p, q, r) of initial size (2, 4, 3) with parameters and , implying

Fig. 4

Growth of the aggregate. Simulation of an aggregate (p, q, r) of initial size (2, 4, 3) with parameters and , implying

Plotting , , and (not shown) for the solutions depicted in Fig. 4 suggests , as . This is in agreement with the following formal result.

Proposition 1

With the notation of Theorem 2, if , then there exists a formal approximation of a solution of (12) of the formwithThe approximation is (from a formal point of view) unique, including the choice of the exponents of t, among solutions with polynomially or exponentially growing aggregate size r.

Remark

The precise meaning of the uniqueness statement is best understood from the proof below. The attempt to find an exponentially growing approximation for a solution by formal asymptotic expansions fails. A polynomial ansatz, on the other hand, leads to the result of the theorem, and the exponents and coefficients are determined uniquely by the asymptotics.

Another result of the parallel work (Delacour et al. 2020) is a rigorous version of Proposition 1. It confirms the existence of solutions with the asymptotic behavior stated in the theorem, and it shows that they are locally attracting. The proof involves advanced techniques from dynamical systems theory such as Poincaré compactification and geometric singular perturbation theory.

In the formal arguments below we shall assume that the solution components behave asymptotically as either like with or like with . This includes the assumption that solution derivatives also behave like the derivatives of these functions.

Since holds for admissible states, when r(t) tends to infinity, then also q(t) tends to infinity at the same rate, which we write with the sharp order symbol asWith , we write the equations for r and for asSince the right hand sides have to be asymptotically nonnegative by the growth of q and r, taking (30) into account, the first equation implies , and the second implies , i.e.If the growth were exponential, i.e. , , then the exponential growth on the left hand sides of (31) would have to be balanced by terms of the same order of magnitude on the right hand sides, i.e. . Then the quadratic negative term on the right hand side of the first equation in (12) could not be balanced by any of the linear positive terms, and would eventually drive p to negative values. This contradiction excludes exponential growth.

For polynomial growth, i.e. , the term in (31) needs to be balanced, implying . In the equation for q in (12), are small compared to . Therefore it is necessary that , implying and, thus, . This justifies the ansatz (28) with the addition . Substitution into the differential equations and comparison of the leading-order terms gives equations for the coefficients:This system can be solved explicitly by first noting that the first two equations imply and, thus, . Using this in the third and fourth equation gives a linear relation between and . This again can be used in the fourth equation to write as a linear function of . The division of the first equation by then gives the formula for in (29). The positivity of is a consequence of (27).

For all the results so far the positivity of the rate constant for breaking cross-links has been essential. Therefore it seems interesting to consider the special case separately. It turns out that the dynamics is much simpler. The aggregate size always grows linearly with time.

Theorem 3

Let , , and . Let satisfy (14). Then the solution of (12) satisfies

For the right hand sides in (12) depend only on p and , meaning that these two variables solve a closed system:The unique nontrivial steady state can be computed explicitly. We prove that it is globally attracting by constructing a Lyapunov functional. Let andFor each point there is a unique value of such that . Therefore the Lyapunov functionis well defined and definite in the sense with equality only for . It remains to prove that the flow on is strictly inwards. For example, for the left boundary part,where the first inequality follows from , and the second inequality from and from the fact that the last parenthesis vanishes for , since is an equilibrium. Analogously, for the right boundary part,for the top boundary part,and for the bottom boundary part,This shows that the Lyapunov functional is strictly decreasing along solutions, completing the proof of global asymptotic stability of the equilibrium.

The linear growth of q and r follows from

This result shows that the breakage of cross-links has somewhat contradictory effects, depending on the parameter regime: It can speed-up the aggregation dynamics, producing a quadratic rather than linear growth of the aggregate size (Case 3 of Conjecture 1 and Proposition 1). This is linked to the fact that it allows the aggregates to rearrange in a more compact way. On the other hand, it may slow down the dynamics, such that the aggregate only reaches a finite size (Case 1) or even disintegrates completely (Case 2), again compared to the linear growth without cross-link breaking (Theorem 3).

Comparison with experimental data: discussion

Comparison with experimental data

There are only limited options for a serious comparison of the theoretical results with experimental data. We shall use the data shown in Fig. 5, which have been published in Zaffagnini et al. (2018). It provides observed numbers of aggregates in dependence of ubiquitin for a fixed concentration of p62. Our results do not permit a direct comparison with this curve, which would require modelling of the process of nucleation of aggregates. However, the data provide at least some information about concentration levels of ubiquitin and p62, such that stable aggregates exist.

Fig. 5

Number of aggregates in terms of [Ubi] (or more precisely ) at fixed M (Zaffagnini et al. 2018). Average and SD among three independent replicates are shown. The dashed line represents a fitted sigmoidal (more precisely, logistic) function, centered around M. Note that here p62 monomers are counted. Under the assumption that p62 only occurs in oligomers of size n we have [p62] = n[p62]. The regression coefficient R measures the quality of the fit

For meaningful quantitative comparisons with these scarce data we need to reduce the number of parameters in our model. As a first step, we fix the value of the size of p62 oligomers, following Zaffagnini et al. (2018) where values between 5 and 6 for GFP-p62 have been found (although we note that in Wurzer et al. (2015) an average of about has been reported for mCherry-p62 in vitro). Since the experiment in Zaffagnini et al. (2018) has been carried out with a fixed p62 concentration M, this corresponds to an oligomer concentration of M.

Concerning the rate constants, we make the assumption that the binding and, respectively, the unbinding rate constants are equal, i.e. and . This will allow to express all our results in terms of one dissociation constant .

From Fig. 5 we conclude that for an oligomer concentration of M the growth of stable aggregates requires a minimal cross-linker concentration [Ubi] roughly between 0.6 and M (M, where M and M are rough estimates of the inflection point and, respectively, the spread of the curve). According to the results of the preceding section, these values should correspond to situations with either or , depending on the question, if the equilibrium aggregate sizes of Case 1 in Conjecture 1 are large enough to be detected in the experiment, or if we need to be in Case 3 of growing aggregates. Therefore, with the above assumptions, with , , and with (26), (27), we obtain for :and for :Solving these equations for with , M, and with [Ubi] between M and M, gives estimates for between M and M for , and between M and M for . So we claim that at least the order of magnitude is significant. It differs by three orders of magnitude from published data on the reaction between ubiquitin and the UBA domain of p62 [M (Long et al. 2010)]. This should not be so surprising, since in the context of growing aggregates the reactions can be strongly influenced by avidity effects.

Discussion

We return to Conjecture 1, where the long-time behaviour is described in terms of the value of the parameter defined in (20). With the simplifying assumptions on the reaction rate constants from above, the statements of the conjecture are depicted in Fig. 6 for the fixed values and M (motivated by the estimates above) in a bifurcation diagram in terms of the concentrations [Ubi] and . In the derivation of the model we have stated the assumption that [Ubi] and are of the same order of magnitude, meaning that their ratio has been kept fixed when taking asymptotic limits (). The bifurcation diagram gives some indication of what happens when they take extreme values. The right hand sides of (33) and (34) indicate that smallness of [Ubi] or can be compensated by largeness of the other, to obtain stable or growing aggregates. However, the occurrence of [Ubi] on the right hand sides produces an unsymmetry. If aggregates are unstable, no matter how large [Ubi] is. Similarly, polynomial growth never happens for .

Fig. 6

Bifurcation diagram corresponding to Conjecture 1 for , M with the bifurcation curves given by (33) and (34)

There is a significant uncertainty concerning the oligomer size n, which has so far been assumed to be 5, according to observations in Zaffagnini et al. (2018). Actually, a distribution of oligomer sizes should be expected in the experiments of Fig. 5 with the occurrence of much larger oligomers. An extension of our model to a situation with a given size distribution of free oligomers would be feasible. Because of the fact that the rate of adding oligomers to the aggregate is independent of the oligomer size [see (3)] it would then be reasonable to assume that also within the aggregate the relative oligomer size distribution is the same as for free oligomers. Without having carried out the computations in detail, we expect a model of the same form as ours, where quantities depending on the oligomer size have to be replaced by suitably computed mean values. There does not seem to be any reason for a different qualitative behavior of the modified model. As an indication of the effects of varying the oligomer size the computation of from (33) has been repeated for a range of values of n between and . The results are depicted in Fig. 7, which shows that the predicted values of might be larger by up to an order of magnitude compared to the case , but still small compared to Long et al. (2010), if larger oligomer sizes are considered and is relevant. The asymptotic behaviour for large oligomer sizes is easily seen to be . On the other hand, if is relevant, the value of becomes smaller by up to an order of magnitude for large oligomers with the asymptotic behaviour . These asymptotic behaviors are easily seen from the explicit solutions of the quadratic equations (33) and (34).

Fig. 7

The dissociation constant determined from (33) (solid line) and (34) (dashed line), depending on the p62 oligomer size n. Ubiquitin and p62 oligomer concentrations from Fig. 5 at the onset of aggregation: M, M

Conclusion

In this article, we have proposed an ODE model for the growth and decay of aggregates of p62 oligomers cross-linked by ubiquitin chains. Under the assumption of unlimited supply of free oligomers and cross-linkers we found three possible asymptotic regimes: complete degradation of aggregates, convergence towards a finite aggregate size, and unlimited growth (quadratic in time) of the aggregate size. In the latter case, growing aggregates are asymptotically tightly packed with the maximum number of cross-links. These statements are supported by a mixture of explicit steady state computations, formal asymptotic analysis, and numerical simulations. The three regimes, which can be separated explicitly in terms of the reaction constants, have been illustrated by the simulation results. Rigorous proofs of some results on the long-time behaviour in the three regimes are the subject of the parallel work (Delacour et al. 2020).

A comparison of the theoretical results with data from Zaffagnini et al. (2018) has provided an estimate for the dissociation constant of the elementary reaction between ubiquitin and the UBA domain of p62 in the context of growing aggregates.

There are several possible extensions of this work. A limitation of the original discrete model is that the description of aggregates by triplets (i, j, k) is very incomplete. Typically, very different configurations are described by the same triplet. For example, we could imagine very homogeneous or very heterogeneous aggregates, i.e. fully packed in certain regions and very loose in others. Reaction rates will strongly depend on the configuration, including information about the geometry of the aggregate. In principle one can imagine an attempt to overcome these difficulties based on a random graph model (Frieze and Karonski 2015), but the resulting model describing probability distributions on the sets of all possible aggregate shapes would be prohibitively complex. An intermediate solution would be a more serious approach to finding formulas for quantities like the probability of losing an oligomer, when a cross-link breaks, based on typical probability distributions.

The model (12) describes an intermediate stage of the aggregation process. On the one hand, the large aggregate assumption means that we are dealing with the growth of already developed aggregates, neglecting the nucleation process, which is important for the number of established aggregates. A model of the nucleation process would be based on the discrete representation and it would have to be stochastic. On the other hand, we neglect two effects important for a later stage of the process. The first and obvious one is the limited availability of free p62 oligomers and ubiquitin cross-linkers. It would be rather straightforward to incorporate this into the model, however at the expense of increased complexity. It would also eliminate the dichotomy between the Cases 1 and 3 of Conjecture 1 since unbounded growth would be impossible. For relatively large initial concentrations of free particles, one could imagine a two-time-scale behaviour with an initial quadratic growth and saturation on a longer time scale. The other effect, which is neglected here but definitely present in experiments, is coagulation of aggregates. This is the subject of ongoing work, based on the PDE model (38) derived in the Appendix and enriched by an account of the coagulation process. A more complete model for a distribution of aggregates, including nucleation and coagulation, would permit a more substantial comparison with experiments as in Zaffagnini et al. (2018), whose output should be extended to provide information about aggregate size distributions. A new experimental challenge suggested by our model is to detect the oligomer-to-crosslinker ratio in aggregates, which could be compared to our theoretical predictions.

More generally, it is very common that protein aggregates contain different types of molecules. In particular, the occurrence of ubiquitinated aggregates has been observed in pathological conditions [e.g. in neurodegenerative diseases (Donaldson et al. 2003) or in diabetes (Kaniuk et al. 2007)]. Existing mathematical models of these aggregation processes (see, e.g. Prigent et al. 2014) could be extended by the approach of this work. Finally, various cytoskeletal structures can be seen as heterogeneous protein aggregates, for example actin filament networks with filamin, fascin, or Arp2/3 complex as cross-linkers (Vinzenz et al. 2012), or focal adhesions, containing integrins, talin, vinculin, and several other protein species (Zamir and Geiger 2001). The latter seems particularly interesting because of the dependence of the growth behaviour on mechanical stimuli (Riveline et al. 2001). Again the approach of this work might provide an alternative basis for mathematical modeling.