Original Research ARTICLE
A Model of Autophagy Size Selectivity by Receptor Clustering on Peroxisomes
- 1Department of Physics and Atmospheric Science, Dalhousie University, Halifax, NS, Canada
- 2Department of Physics, Simon Fraser University, Burnaby, BC, Canada
Selective autophagy must not only select the correct type of organelle, but also must discriminate between individual organelles of the same kind so that some but not all of the organelles are removed. We propose that physical clustering of autophagy receptor proteins on the organelle surface can provide an appropriate all-or-none signal for organelle degradation. We explore this proposal using a computational model restricted to peroxisomes and the relatively well characterized pexophagy receptor proteins NBR1 and p62. We find that larger peroxisomes nucleate NBR1 clusters first and lose them last through competitive coarsening. This results in significant size-selectivity that favors large peroxisomes, and can explain the increased catalase signal that results from siRNA inhibition of p62. Excess ubiquitin, resulting from damaged organelles, suppresses size-selectivity but not cluster formation. Our proposed selectivity mechanism thus allows all damaged organelles to be degraded, while otherwise selecting only a portion of organelles for degradation.
Macroautophagy (hereafter autophagy) can degrade large subcellular substrates such as organelles and pathogenic bacteria [1–3]. Selective autophagy can target specific damaged or surplus substrates , and defects in this process often lead to human disease . To degrade targeted substrates, the downstream stages of autophagy include formation of and recruitment to phagophores, and formation and maturation of the autophagosome. The initial targeting of substrates involves an “eat-me” signal such as ubiquitin and subsequent recruitment of autophagy receptor proteins, with only some receptor types recruited for a given type of organelle .
Selective autophagy involves an all-or-none response, where each substrate is either targeted or not for degradation, leading to a relatively very high or relatively low respective rate of degradation by the autophagy system. Such all-or-none responses can be produced through cooperative effects [6, 7]. Clustering is a cooperative mechanism that can generate a qualitative all-or-none response, as illustrated by the regulation of bacteriophage lysis timing where lysis follows only after the collective formation of holin-clusters on the cell surface [8, 9]. The involvement of protein clusters in autophagy is supported by reports of domains of distinct receptor proteins on bacteria targeted for xenophagy [10–12].
Autophagy receptor proteins are relatively well characterized for mammalian peroxisomes. Peroxisomes are essential and dynamic cellular organelles with functions including metabolism of hydrogen peroxide and oxidation of fatty acids . Peroxisomes are approximately spherical with diameters ranging from ~0.1−1μm . There can be hundreds of peroxisomes in a single mammalian cell . In both mammals and yeast significant peroxisome degradation is through autophagy, known as pexophagy [16, 17]. In mammalian autophagy, peroxisomes are targeted for degradation by exogenous ubiquitin labeling . However, the receptor protein NBR1 is also both necessary and sufficient for pexophagy  while the receptor p62 significantly contributes though is not essential [18, 19]. The J region of NBR1  allows NBR1 to anchor directly to organelle membranes.
We propose that a sufficient all-or-none signal for autophagy selectivity can be provided by the initial formation of receptor clusters on substrate surfaces. Specifically, for mammalian peroxisomes, we hypothesize that NBR1 can form lateral clusters once it is associated with membranes, and that an NBR1 cluster on a peroxisome is necessary for downstream degradation while the absence of an NBR1 cluster prevents downstream degradation. In this paper we explore this cluster-selectivity hypothesis in the context of mammalian peroxisomes. Small clusters of receptor proteins, if initially placed on every subcellular organelle, will subsequently grow and shrink due to receptor exchange between organelles . With such initial cluster placement, selectivity would only emerge at late times when only a few organelles are left with clusters of receptors. We consider here the selective initial formation of NBR1 clusters, the role of ubiquitin on NBR1 recruitment, and the possible effect of p62 on NBR1 cluster formation. We take a computational modeling approach to investigating our hypotheses, so that we can quantitatively test the self-consistency of our model .
Our computational model is deliberately simple with respect to the complex biological regulation of cellular processes. Our goal is simply to demonstrate a physically viable mechanism by which autophagy selectivity might operate, consistent with the known biophysical functions of the peroxisomal receptor proteins NBR1 and p62. With this computational approach, we both pose and address three essential questions about autophagy selectivity: (1) what mechanism can provide an all-or-none signal to target individual organelles for degradation by autophagy, (2) how can this mechanism respond flexibly to organelle damage, and (3) what needs to be measured in the lab to better characterize mechanisms of autophagy selectivity? We also identify robust qualitative implications of our hypotheses that can be confronted with experiments, even in the face of considerable parameter uncertainty [22, 23].
2.1. NBR1 Model Description
We start with an NBR1-only model, since NBR1 is the only peroxisomal receptor that is both necessary and sufficient for mammalian pexophagy . The results of our NBR1 model development are presented here without mathematical details (see details below in Computational Methods). Our model for NBR1 dynamics has three aspects: NBR1 association with the peroxisome surface, formation of NBR1 clusters from freely-diffusing NBR1 on the peroxisome surface, and the growth or shrinkage of NBR1 clusters once they have formed. These aspects are illustrated in Figures 1A–C.
Figure 1. Model processes. (A) NBR1 (red circles) in the cytosol surrounding a peroxisome (large purple sphere) associates with the peroxisome membrane via transient association with ubiquitin (blue squares) bound to the peroxisome membrane (association rate in Equation 2). NBR1 on the peroxisome membrane can then dissociate from the peroxisome membrane (dissociation rate in Equation 4). (B) With sufficient NBR1 present on the peroxisome membrane, a cluster of NBR1 nucleates (required concentration in Equation 10). The number of NBR1 in a cluster can increase or decrease following nucleation (Equation 1). Peroxisomes colored green harbor an NBR1 cluster and may be selected for degradation by autophagy. (C) Ostwald ripening results in the growth of larger clusters, and shrinking and evaporation of smaller clusters (Equation 1 describes cluster size dynamics). (D) Under “disuse” conditions that cause the number of ubiquitin on all peroxisomes to increase, more NBR1 will associate with peroxisome membrane, facilitating cluster nucleation, and the selection of larger peroxisomes for degradation by autophagy. (E) Damage may cause the number of ubiquitin on a subset of peroxisomes to increase, leading to increased NBR1 association on those peroxisomes, possible cluster nucleation, and selection of those damaged peroxisomes for degradation by autophagy.
Peroxisomes have surface-displayed ubiquitin [18, 24], and NBR1 have a ubiquitin binding region (UBA) [25–27], so NBR1 can associate with peroxisomes by attaching to the surface-displayed ubiquitin. However, the dissociation constant of the UBA region is much larger than estimated NBR1 concentrations—strongly suggesting that NBR1 will only transiently associate with ubiquitin (see details in Computational Methods). This leaves the amphipathic J domain of NBR1  as the dominant mode of lasting association, following much lower dissociation constants for amphipathic helices.
We model NBR1 recruitment to the peroxisome surface as diffusion-limited arrival to and transient association with the ubiquitin on the peroxisome surface via the UBA domain of NBR1, immediately followed by long-lived association with the peroxisome membrane via the J region of NBR1. This is quantitatively described by a standard equation for diffusion-limited association to absorbing targets on a sphere, see Equation (2) below.
Once associated with the peroxisomal membrane, our model allows NBR1 to self-associate into homo-oligomeric clusters, motivated by the coiled-coil domains of NBR1 [19, 28], the importance of NBR1 to aggregation formation , and generic aggregation phenomenon driven by non-specific interactions [30–33]. Formation of clusters from many freely-diffusing individual NBR1 already on the peroxisome surface occurs at a critical concentration of individual NBR1. To determine this critical concentration for cluster formation we require sufficient NBR1 to both form a cluster (condition 1) and leave behind enough individual NBR1 to prevent the cluster from immediately evaporating (condition 2). The critical concentration is the lowest concentration at which both of these conditions can be satisfied, given by Equation (9) below, which results in the formation of a cluster.
Our model of NBR1 cluster size increase and decrease is in line with the standard physical picture for such processes, known as Ostwald ripening [34, 35]. Clusters described by Ostwald ripening will shrink if they are below some threshold size, and grow if they are above the threshold size, with this threshold size growing in time [34, 35]. The continuous evolution of an individual cluster size with time t is given by
which is derived below leading up to Equation (15). Nclust is the number of NBR1 in a cluster, and R is the radius of the peroxisome harboring the cluster. The total cellular NBR1 which is not associated with peroxisomes is in a shared cellular pool, and it is through this pool that NBR1 can exchange between different peroxisomes. Two variables, a(t) and w(t), are time-dependent combinations of variables such as the cellular concentration of NBR1, or the number of freely-diffusing NBR1 on the peroxisome under consideration. The remaining parameters (σ∞, ν, and b) are constants, independent of time or the peroxisome under consideration.
2.2. NBR1 Results
Using our quantitative NBR1 model, we examine NBR1 dynamics on individual peroxisomes. This includes the behavior of the number of associated NBR1 with time, the timing of the formation of NBR1 clusters at a critical concentration, and the growth of clusters after they have formed.
2.2.1. NBR1 Clusters form Promptly and Allow Much More NBR1 to Associate with Peroxisomes
When peroxisomes are placed in a medium with a bulk concentration of NBR1 (at time t = 0), Figure 2A shows how the recruitment of NBR1 by ubiquitin leads to an approximately linear increase of freely diffusing surface-associated NBR1 (dashed blue line). At very early times there are no NBR1 clusters (dotted red line), so the total NBR1 on the peroxisome (solid black line) consists only of NBR1 freely diffusing on the surface. However, once the surface NBR1 concentration surpasses the critical concentration of nucleation, at approximately t = 10 s, a cluster is nucleated (top process in Figure 1B) and very quickly grows in size. At late times almost the entire population of surface-associated NBR1 is in clusters, while the population of freely-diffusing NBR1 decreases somewhat after cluster nucleation. This decrease is due to the reduced supersaturation necessary to maintain a larger cluster as compared to the supersaturation necessary to nucleate a smaller cluster (see Equation (7) below and ).
Figure 2. (A) The amount of NBR1 associated with a peroxisome of radius R = 0.25 μm vs. time. The initial (t = 0) bulk concentration is [NBR1] = 10 μm−3, with no surface NBR1. At approximately time t = 10 s, a cluster is nucleated and grows thereafter. Cluster growth (dotted red line) allows the total NBR1 (solid black line) to greatly exceed the saturation level exhibited by the free NBR1 (dashed blue line). (B) The total amount of NBR1 vs. time is shown for each of 10 peroxisomes, with radii exponentially distributed between R = 0.05 μm and R = 0.50 μm. The largest peroxisomes have the most NBR1, and the smallest the least NBR1. For each peroxisome cluster nucleation is indicated by a circle, evaporation by a square, and a thicker line indicates when a cluster is present. Only the largest five peroxisomes form clusters, starting around t = 20 s. For all but the largest peroxisome, the clusters nucleate, grow, and then shrink again as the clusters coarsen . Inset of (B) shows the NBR1 surface concentration on peroxisomes vs. time—the concentrations for the different peroxisomes, all of different radii, are equal except when a cluster is present. (C) The size of the smallest peroxisome with a cluster, Rclust, vs. time shows how at early and late times only the largest peroxisomes have clusters. At intermediate times, clusters form and remain on all but the smallest peroxisomes. Rclust reaches a minimum at intermediate times, indicating the radius of the smallest peroxisome to be occupied by a cluster at any time, R*. The legend indicates the initial bulk NBR1 concentration—higher concentrations lead to lower R*.
The behavior of freely-diffusing and cluster-bound NBR1 illustrates how the existence of clusters can lead to a prompt all-or-none response—either a cluster comprised of a large number of NBR1 is present, or a cluster is not present and only a much smaller number of NBR1 are present on the peroxisome surface. As a result of cluster formation, the total NBR1 on peroxisomes can greatly exceed the maximal (saturated) level of freely-diffusing NBR1. After only 1,000 s, in Figure 2A, we observe more than an order of magnitude excess of NBR1 in the cluster (dotted red line) than in the highest surface concentration (dashed blue line). NBR1 cluster formation provides a prompt all-or-none response that can be a mechanism for pexophagy selectivity.
2.2.2. Largest Peroxisomes Acquire NBR1 Clusters First, and Lose Them Last
The equations describing our quantitative model of NBR1 dynamics all depend on the peroxisome radius R—e.g., the rate at which NBR1 associates with the peroxisome surface; the critical concentration for cluster formation; and the change in cluster size in time. [See Equations (2), (9), and (15) below, respectively.] Because the NBR1 dynamics depend quantitatively on the peroxisome size, the behavior of NBR1 on a particular peroxisome, both freely-diffusing and cluster-bound, is affected by the radius of that particular peroxisome. They are also affected by the radii of other peroxisomes, since the NBR1 populations on all peroxisomes share the same cellular pool of NBR1.
In Figure 2B, there is a system of 10 peroxisomes, with polydisperse radii between R = 0.05μm and R = 0.50μm. (Larger numbers of peroxisomes behave similarly, but are not as easily visualized. Each line in the figure corresponds to the number of NBR1 on a different individual peroxisome). Initially, the freely-diffusing NBR1 on all 10 peroxisomes increases approximately linearly, similar to the behavior of Figure 2A. Larger peroxisomes are able to recruit more NBR1, and so in Figure 2B the largest peroxisome has the most NBR1, followed by the next largest peroxisome, and so on. Once the largest peroxisome reaches the critical concentration of freely-diffusing NBR1, soon after t = 10 s, it nucleates a cluster. This is followed by the next largest peroxisome nucleating a cluster, etc, until the five largest peroxisomes have a cluster, after which no further cluster nucleation occurs. The NBR1 pool not associated with peroxisomes is depleted as the nucleated NBR1 clusters grow and sequester NBR1; this reduction prevents smaller peroxisomes from reaching the critical NBR1 concentration and forming clusters.
Figure 2B shows that large peroxisomes form clusters, while small peroxisomes do not; and that clusters on larger peroxisomes grow earlier and shrink later than clusters on smaller peroxisomes. This implies that peroxisomes below a certain size will not have clusters.
The clusters which have formed on the five largest peroxisomes then proceed to compete for material—the smallest cluster, which is on the smallest peroxisome with a cluster, shrinks until it evaporates (bottom process in Figure 1B), followed by the next smallest cluster, until there is only a single cluster remaining—this is shown schematically in Figure 1C. This competition between clusters for NBR1, mediated by bulk diffusion between clusters, is known as Ostwald ripening [34, 35]. Only the larger drops retain clusters at late times .
The initial nucleation of clusters is prompt and selects for larger peroxisomes. Furthermore, the initial nucleation of clusters is orders of magnitude faster than the slow resolution of clusters through Ostwald ripening dynamics.
2.2.3. Long Lifetime for NBR1 Clusters on Even the Smallest Occupied Peroxisomes
Figure 2C examines how the smallest peroxisomal radius that has a cluster, Rclust, changes in time. We consider an exponentially distributed polydisperse distribution of peroxisome radii between R = 0.05 μm and R = 0.50 μm, using a system of 100 peroxisomes for good size resolution. We track the radius of the smallest peroxisome that has a cluster, Rclust, vs. time. Changes in Rclust are due to the combined effect of nucleation of new clusters and evaporation of existing clusters. Rclust shows that large peroxisomes are selected for cluster formation and growth.
At around t = 10 s cluster nucleation begins with the largest peroxisomes, followed by nucleation on progressively smaller peroxisomes. Eventually cluster formation halts, even though not all peroxisomes are occupied by a cluster. As the larger peroxisomes formed clusters first, the end of cluster formation defines the radius of the smallest peroxisome occupied by a cluster at any time, R*. Peroxisomes of radius R* and larger harbor a cluster at some time, while smaller peroxisomes are always unoccupied. Later, around t = 600 s, clusters begin to evaporate from the smaller occupied peroxisomes—and this “Ostwald ripening” causes Rclust to slowly increase with time.
Each curve in Figure 2C represents a different NBR1 concentration, as indicated in the legend. Higher NBR1 concentrations lead to earlier nucleation and evaporation, a lower Rclust at all times, and a lower R*. In every case, we see that there is more than a decade in time where the smallest occupied peroxisome retains a cluster. During this period the size of that cluster is first growing and then shrinking. Clusters on larger peroxisomes than R* survive even longer, as illustrated in Figure 2B.
2.3. Size Selectivity and Average Peroxisome Size
Under our hypothesis that NBR1 clusters are required for pexophagy selectivity, the absence of NBR1 clusters on the smallest peroxisomes suggests some size-selectivity (determined by R*). While our model does not directly address downstream pexophagy processes, if they are either fast enough to respond immediately to cluster formation (which progresses from the largest to the smallest peroxisomes) or slow enough to allow for evaporation of smaller NBR1 clusters (which progresses from the smallest to the largest peroxisomes) then size-selectivity could be even more substantial.
We can explore some downstream effects of size-selectivity, using R*. We explore three simplified cases: maximal, partial, or zero size selectivity. For maximum size selectivity all of the largest peroxisomes, defined as those with a radius greater than R*, are degraded. Since larger peroxisomes are expected to harbor larger clusters (Figure 2B), maximum size selectivity corresponds to the hypothesis that downstream processes first select peroxisomes with the largest NBR1 clusters. For partial size selectivity, a fraction of all peroxisomes with radii greater than R* are degraded. This corresponds to the hypothesis that downstream processes randomly select peroxisomes with NBR1 clusters, but do not otherwise differentiate between them. For zero size selectivity, peroxisomes are selected randomly for degradation. This corresponds to the hypothesis that NBR1-clusters play no role in autophagy selectivity.
A common measure of pexophagy activity is relative changes in the amount of catalase [18, 19, 37], since catalase is an abundant peroxisomal marker. We will show that size-selectivity can significantly affect a volumetric measure such as total catalase fluorescence intensity, even when the number of degraded organelles remains unchanged.
2.3.1. Size Selectivity Significantly Affects Catalase Abundance
Maximum, partial, and zero size selectivity are illustrated in Figure 3A for exponentially distributed peroxisomal radii between Rmin = 0.05 μm and Rmax = 0.5 μm. The red shaded region corresponds to a fixed fraction of degraded peroxisomes under the three scenarios of maximal, partial, or zero selectivity (as indicated). For maximum size selectivity, all peroxisomes with radii greater than R* are degraded. For partial size selectivity, only peroxisomes with radii greater than R* are degraded, however only a fraction p of the peroxisomes in this size range are randomly degraded. For zero size selectivity, peroxisomes are selected for degradation randomly, so that peroxisomes of different size are selected in proportion to their abundance. Maximum and zero selectivity correspond to limiting cases of partial selectivity, with p = 1 or , respectively.
Figure 3. (A) Maximum size selectivity in pexophagy is achieved when all of the largest peroxisomes, of radius R > R*, are degraded (shaded red area) from a distribution of peroxisome sizes (solid black line, illustrating an exponential distribution of abundance P(R) vs. radius R). Partial size selectivity is achieved when a fraction p of the larger peroxisomes, with R > R*, are degraded. Zero size selectivity is achieved when peroxisomes are randomly chosen for degradation, independent of size. (B) The remaining peroxisomal volume fraction vs. the remaining number fraction r, after various amounts of zero selectivity degradation (dashed diagonal red line) or maximal selectivity degradation (solid black line, given by Equation 16). Partial selectivity is equivalent to maximum selectivity a fraction p of the time, and so lies on a line interpolating between maximum selectivity with the same R* and no autophagy (at r = 1). One such interpolating line (dotted blue) is illustrated.
We consider a scenario where a certain fraction of the total number of peroxisomes are degraded, with a fraction r remaining. For the same number of peroxisomes removed, maximum, partial, and zero size selectivity will each result in a different fraction of the total peroxisomal volume being removed upon degradation. We calculate the volume fraction remaining as peroxisomes are degraded (shown in Computational Methods), and the results of these calculations are shown in Figure 3B. For zero size selectivity, the remaining peroxisomal volume is strictly proportional to the remaining peroxisomal number (dashed red line). For maximum selectivity, the remaining volume fraction decreases sharply as the number fraction is decreased (solid black line). This is because only the largest peroxisomes are targeted when size-selectivity is maximal. Partial selectivity (dotted blue line) interpolates between maximum selectivity and no autophagy (at r = 1). Only one interpolating line is shown for illustrative purposes.
We see that the reduction in peroxisome volume, and the corresponding reduction in catalase intensity, is a combination of the number of peroxisomes removed and the size selection of those peroxisomes. Many more peroxisomes need to be removed with zero selectivity to match the same volume reduction at maximum selectivity. Since a change of catalase intensity is a proxy for a change in peroxisomal volume, our results indicate that catalase intensity can be affected by either changing the number of peroxisomes degraded or by changing the size-selectivity.
2.4. p62 Model Description
Experimental p62 inhibition with siRNA causes an increase in the catalase signal . This relative increase of catalase has been previously interpreted as a decrease in the pexophagy rate due to p62 inhibition , but any change in size-selectivity could confound this interpretation. To explore this possibility, we now extend our model to consider how p62 could affect size selectivity of peroxisomes. The results of our p62 model development are presented here without mathematical details (see details below in Computational Methods). Our model for p62 dynamics has two aspects: p62 association with and recruitment to peroxisomes, and p62 inhibition of NBR1 clusters.
Similar to NBR1, p62 can in principle bind to ubiquitin associated with peroxisome membranes using its UBA domain [38, 39]. However (details in Computational Methods), the dissociation constant is too large to expect significant p62 binding to ubiquitin. Nevertheless, the p62 PB1 domain has a strong affinity to other PB1 domains , and so p62 can bind to membrane-associated NBR1 through its PB1 domain. Thereafter, associated p62 can form filaments through PB1-PB1 interactions [41, 42].
p62 recruitment to NBR1 on the peroxisome surface is modeled as diffusion-limited arrival to and association with NBR1 on the peroxisome membrane, or to p62 already associated with NBR1. These are both quantitatively described by a standard equation for diffusion-limited association to absorbing targets on a sphere (Equation 3 below).
Polymer physics leads us to hypothesize that polymeric chains of p62 associated with NBR1 could reduce NBR1 self-association through steric repulsion. Such a repulsion arises as the entropic contribution to the free-energy of the polymers is decreased when brought close together . Such steric repulsion of membrane associated proteins can prevent growth of protein clusters , can lead to cluster segregation , and can even inhibit phase separation of associated lipids . To explore the consequences of a strong steric repulsion between p62-associated NBR1, our model does not allow NBR1 that is associated with p62 to participate in cluster formation or growth.
2.5. p62 Results
In our quantitative model a p62 filament on a membrane-associated NBR1 prevents that NBR1 from participating in cluster formation or growth—which should reduce NBR1 numbers. Despite this cluster-inhibition mechanism, NBR1 and p62 both contain a LIR domain, which interacts with the machinery for the formation of autophagosomes [25, 47]. Given these opposing effects, how does the number of membrane-associated LIR domains change as the p62 concentration is varied? How does greater or lesser amounts of surface ubiquitination affect the function of p62?
2.5.1. p62 Inhibits NBR1 Clusters
For an individual peroxisome, we show in Figure 4A how the total steady-state number of LIR on a peroxisome changes as both the global p62 concentration and the peroxisomal number of ubiquitin are varied. For larger p62 concentration (the purple/red region) NBR1 clusters are suppressed, dramatically reducing LIR content. Within this non-clustering regime the LIR content slowly increases with increasing p62. For smaller p62 concentration (the yellow region) we see a striking increase in the LIR count due to formation of an NBR1 cluster. The clustering regime has approximately 10 × the number of LIR domains as the non-clustering regime. Within the clustering regime, increasing p62 decreases the LIR count by decreasing the cluster size.
Figure 4. (A) The steady-state number of LIR domains on a single peroxisome after equilibration, given by the sum of the number of NBR1 and p62, as cellular [p62] (in μm−3) and number of ubiquitin are varied. R = 0.25 μm, [NBR1] = 10 μm−3, and V = 100 μm3. NBR1 clusters are present below the dashed line. Increasing ubiquitin increases the LIR count, and so does increasing [p62] when no clusters are present. However, increasing [p62] suppresses cluster formation and decreases the LIR count within the clustering regime. (B) The radius of the smallest peroxisome with a cluster, Rclust, vs. time, following induction of cluster formation after increase of the ubiquitin coefficient at t = tstep. [p62] is varied as indicated, and we otherwise use the same conditions as in Figure 2C. Initially the largest peroxisomes form clusters, then for a decade of time no more clusters are formed, followed by cluster evaporation from the smallest occupied peroxisomes. Increasing [p62] leads to more selective cluster formation, with only the larger peroxisomes occupied by clusters. (C) Lines indicate the boundary between regimes with at least one NBR1 cluster on some peroxisome (above) and those with no clustering on any peroxisome (below), as both cellular [p62] and ubiquitin are varied as indicated. Different lines correspond to cellular [NBR1] as indicated by the legend. Other conditions are the same as in Figure 2C. For higher ubiquitin, lower [p62], or higher [NBR1] we observe clusters. Even for higher [p62], sufficiently high ubiquitin will still lead to NBR1 clusters—though note the logarithmic scale. (D) After the ubiquitin coefficient n0 is suddenly raised (as indicated by the legend), the data shows the radius of the smallest peroxisome with a cluster at any time, R*, vs. cellular p62 concentration. The disuse and damage scenarios correspond to ubiquitin being increased on all peroxisomes (“all,” diamonds and triangles) or increased only on a single peroxisomes (“one," circles and squares), respectively. The damage scenario leads to significantly smaller R*, as does a larger final ubiquitin coefficient. While the data points represent an initial ubiquitin coefficient of n0 = 10, the black lines represent initial ubiquitin coefficient of n0 = 25, and their agreement indicates that the initial ubiquitin coefficient is not significant. The number of ubiquitin on a peroxisome is , with R0 = 0.25 μm.
2.5.2. p62 Inhibition of NBR1 Clusters Enhances Size Selectivity
In Figure 2C we showed how there is a time-dependent threshold peroxisome size Rclust - below this size the peroxisomes do not have a cluster, and above this size peroxisomes harbor a cluster. Here we investigate how p62 can change this threshold peroxisome size.
In Figure 4B we show the radius of the smallest peroxisome with an NBR1 cluster, Rclust, vs. elapsed time after the ubiquitin level is simultaneously increased on all peroxisomes. (We use ubiquitin increase to initiate autophagy here since it recruits both NBR1 and p62.) We use the same radius distribution as Figure 2C. The ubiquitin increase induces NBR1 uptake to peroxisome surfaces. We consider various p62 concentrations, as indicated by the legend. Increasing the p62 concentration causes an increase in the size of the smallest peroxisome with a cluster, and limits NBR1 clusters to larger peroxisomes. If only peroxisomes with clusters are degraded by autophagy, p62 increases the threshold size for peroxisome degradation, and thus enhances size selectivity.
2.5.3. Sufficient Ubiquitin Allows NBR1 Clusters to Form on Larger Peroxisomes Despite p62
In Figure 4C we show that sufficiently high p62 concentrations (or low ubiquitin coverage) can completely inhibit clusters on all peroxisomes. Each curve indicates the minimum ubiquitin coefficient that leads to NBR1 cluster formation on any peroxisome, under the same radius distribution as Figure 2C. Below the line, no NBR1 clusters form, while above the line there is at least one NBR1 cluster. The different lines, as indicated by the legend, correspond to different cellular NBR1 concentrations. We see that increasing NBR1 leads to more peroxisomes with clusters.
While sufficiently high ubiquitination allows NBR1 cluster formation at any given p62 concentration, the clusters still only form on larger peroxisomes. This is shown in Figure 4D, where the radius of the smallest peroxisome with an NBR1 cluster, R*, vs. [p62] is shown. With more p62, the radius of the smallest peroxisome occupied by a cluster increases.
2.5.4. Single Organelle Damage has Weaker Size Selectivity than Multi-Organelle Disuse
Figure 4B shows how increasing the p62 concentration restricts NBR1 clusters to larger peroxisomes. For each curve in Figure 4B, the radius of the smallest peroxisome with a cluster reaches a minimum, R*. This is the radius of the smallest peroxisome that ever has a cluster.
In Figure 4D we see that R* increases as the p62 concentration increases. We consider two values of the initial ubiquitin coefficient (n0 = 10 with colored points as indicated by the legend and n0 = 25 with corresponding thin black lines) and two values of the final ubiquitin coefficient (n0 = 100 with circles and diamonds, or 150 with squares and triangles, as indicated in the legend). We suddenly change the ubiquitin number on peroxisomes from the low initial value, which does not support clustering, to a high final value that does support clustering. The agreement between points and black lines indicates that the initial ubiquitin coefficient does not significantly affect the clustering. However, the final ubiquitin coefficient does affect clustering.
When we compare the minimum peroxisome size with a cluster, R*, between a “disuse” scenario (shown schematically in Figure 1D), where all peroxisomes increase their ubiquitin (green diamonds and orange triangles in Figure 4D), and a “damage” scenario (Figure 1E), where only one peroxisome increases its ubiquitin (red circles and blue squares in Figure 4D) we see that for damage significantly smaller peroxisomes will form NBR1 clusters. Hence the damage scenario decreases size selectivity, allowing more damaged peroxisomes to be selected with moderate amounts of ubiquitin.
Autophagy selectivity amounts to an all-or-none response, where each substrate is either selected or not selected. Selectivity is thought to be mediated by autophagy receptor proteins. Some autophagy receptor proteins can be observed to form microdomains on substrates [2, 10–12]. Our primary hypothesis is that the presence or absence of receptor clusters on individual organelles provides a necessary all-or-none response for autophagy selectivity, for at least some types of organelles.
There can be hundreds of peroxisomes in a single mammalian cell, and pexophagy can result from either damage of individual peroxisomes or from more generic deproliferation . The signaling protein ubiquitin and the autophagy receptor proteins NBR1 and p62 participate in pexophagy [18, 19, 24]. NBR1 is essential for pexophagy; also essential are the specific domains of NBR1 that govern NBR1 association with itself, with plasma membranes, with ubiquitin, and with autophagosomes. We explore our primary hypothesis by quantitatively modeling the dynamics of NBR1 cluster formation on individual peroxisomes.
With sufficient cellular NBR1 and peroxisomal ubiquitin, our model leads to the formation and growth of NBR1 clusters. Clusters are found to form rapidly on the largest peroxisomes, and subsequently on smaller ones. After cluster nucleation, competitive Ostwald ripening progressively removes initial clusters—starting from the smaller peroxisomes. The smallest peroxisomes (below a critical radius R*) never form clusters. Thus, we find that NBR1 clustering is strongly size-selective—NBR1 clusters are found on large peroxisomes and not on small peroxisomes. Nevertheless, higher NBR1 concentrations or peroxisomal ubiquitin levels can reduce size selection by reducing R*.
While p62 is not essential for pexophagy, experiments show that lower p62 levels increase the cellular abundance of peroxisomal catalase—indicating a significant role in pexophagy . We echo our primary hypothesis with a secondary hypothesis that p62 influences pexophagy by affecting NBR1 cluster formation. Specifically, we hypothesize that p62 association with NBR1 and subsequent p62 filament formation will sterically hinder and inhibit NBR1-NBR1 association. By modeling this hypothesis, we find that p62 can significantly affect size-selectivity and that this alone may be sufficient to explain the catalase response of p62 inhibition.
3.2. NBR1 Clusters and Size-Selectivity
Deosaran et al.  introduced the idea that a “critical mass” of autophagy receptor proteins is necessary to target peroxisomes to autophagosomes. Supporting this idea, when NBR1 exceeds a critical concentration on the peroxisomal surface nucleation of an NBR1 cluster is followed by a sudden and localized further increase in LIR numbers, and provides an all-or-none signal on individual peroxisomes.
We found that cluster formation favors larger peroxisomes because of a lower surface concentration threshold for cluster formation. This is seen in Figures 2B,C, with the first cluster forming on the largest peroxisome and then progressing toward smaller peroxisomes. Subsequent cluster evaporation begins on the smallest peroxisome with a cluster (peroxisomal radius R*) and progresses toward larger peroxisomes . NBR1 clustering provides a consistent signal for autophagy to target and degrade large peroxisomes. This size-selectivity of cluster formation and evaporation, in combination with our hypothesis that NBR1 clusters are an essential part of pexophagy selectivity, is consistent with early reports that larger peroxisomes are preferentially degraded by pexophagy  and that the required proteins for pexophagy depend on peroxisome size .
Cluster formation occurs very quickly, with the number of clusters peaking after an increase in the ubiquitin level in ≲30 s in Figure 4B. Following cluster formation, evaporation of clusters is slower, beginning after 102–104 s in Figure 4B. This evaporation is nevertheless faster than reported for pre-existing clusters  because some clusters are nucleated close to the threshold for evaporation. These cluster formation and evaporation timescales suggest that the selection of peroxisomes by NBR1 cluster formation could occur well within the timescale of mammalian or yeast pexophagy, which take days [16, 50] or hours [51, 52], respectively.
While we hypothesize that NBR1 receptor clusters are necessary for autophagy selectivity, they may not be sufficient. Our scenario of partial selectivity, illustrated in Figure 3A, reflects this possibility.
While NBR1 clusters on peroxisomes and their preference for larger peroxisomes appear to be viable selectivity mechanisms according to our quantitative modeling, they remain hypotheses. We do note that organelle-size dependent activity arising from cluster size modifications has been previously proposed in the context of membrane recycling . Moreover, our proposed mechanism may apply for mitochondria in a cellular anti-viral response mediated by MAVS (also known as CARD or IPS) proteins . Before viral infection, MAVS are not activated, and are spread out on the mitochondrial membrane  and are evenly distributed between mitochondria . Upon viral infection MAVS are activated, and subsequently cluster on the mitochondrial membrane [55, 57]. Significantly, some mitochondria have significant MAVS while other mitochondria have little to no MAVS . Limited data suggests that large mitochondria retain MAVS while small mitochondria do not . The MAVS anti-viral phenomenology appears qualitatively similar to our NBR1 model, and so provides some support for our model behavior.
3.3. Pexophagy Selectivity: Disuse vs. Damage
In our model an increase in ubiquitin number acts as an initial signal that can lead to additional NBR1 and p62 accumulation, possibly followed by autophagic degradation. Ubiquitination of peroxisomes has been shown sufficient to induce pexophagy ; ubiquitin is thought to recruit NBR1, the primary autophagy receptor protein for peroxisomes, to the peroxisomes membrane ; and earlier modeling suggests that an increase in ubiquitin could be a natural and self-correcting response of the cell to requiring fewer peroxisomes . Ubiquitin is known to play a role in the routine import of peroxisome matrix proteins , and is part of the quality control system for damaged proteins on peroxisomes [58–60] that is distinct from the well known role of ubiquitin as a signal for the ubiquitin-proteasome system .
In Figure 4D we explored two extreme cases of increases in the ubiquitin level: a global increase, where ubiquitin levels on all peroxisomes increase, and an increase of the ubiquitin level on a single peroxisome. A global increase (shown schematically in Figure 1D) could be due to the removal of peroxisome proliferators  or a change in growth medium , both of which can result in a decrease in peroxisome numbers as superfluous peroxisomes are degraded. An increase on a single peroxisome (shown schematically in Figure 1E) could be due to damage [62, 63].
For a disuse scenario, size-selectivity would lead to the expectation that larger peroxisomes would be preferentially degraded until the decrease in peroxisome numbers had been achieved and ubiquitin levels reduced  so that NBR1 clusters are no longer formed. This is consistent with observations that larger peroxisomes are preferentially degraded when reducing peroxisome numbers .
For the damage scenario, there is little competition for the NBR1 needed to form clusters since only a few peroxisomes have an elevated ubiquitin level. This results in less size-selection (lower R*) than in the disuse scenario, as seen in Figure 4D. Nevertheless, even in the damage scenario larger peroxisomes more easily form NBR1 clusters, suggesting that larger peroxisomes may need less damage to be selected for autophagy.
3.4. Selectivity with p62
The autophagy receptor protein p62 is important for pexophagy . Inhibition of p62 with siRNA increases total peroxisomal catalase . We have demonstrated that when p62 inhibits NBR1 clustering in our model, increased p62 levels cause greater size-selectivity: the size threshold above which peroxisomes have NBR1 clusters is pushed to larger peroxisomes by increased p62 concentrations. This in itself would be sufficient to explain the catalase increase following p62 inhibition. Volumetric measures of autophagy such as catalase will report the combination of number and volume. To assess the number of organelles targeted by autophagy, the number should be directly assessed.
3.5. Experimental Signatures of NBR1 Clusters
(a) The most striking result of our hypothesis, that NBR1 clusters are necessary for downstream degradation by the autophagy system, is significant size-selectivity. Large peroxisomes will be preferentially degraded over small peroxisomes. One way of measuring this effect would be to measure the fluorescence intensity of a tagged peroxisomal protein, such as catalase, together with the degree of colocalization with a protein associated with autophagosomes, such as LC3 . Our model results indicate that peroxisomes with significant colocalization would have a larger average catalase intensity compared to peroxisomes with little or no colocalization.
(b) Our model also indicates that formation of NBR1 clusters will significantly affect the number of NBR1 associated with peroxisomes of similar size. For the parameters of Figure 4A, the differences are approximately 10-fold. While surface ubiquitin concentrations may differ between peroxisomes of similar sizes, and so could determine which peroxisomes have clusters, ubiquitin alone appears unlikely to be able to directly affect non-cluster NBR1 to the same extent .
With NBR1 clustering, we have shown that peroxisomes will either have a cluster and have a large amount of NBR1, or not have a cluster and have a small amount of NBR1. After pexophagy is induced, e.g., by the removal of peroxisome proliferators, the fluorescence of NBR1 colocalizing with catalase should therefore have a bimodal distribution. Qualitatively, the all-or-none colocalization of NBR1 with peroxisomes suggested by, e.g., Figure 5 of Deosaran et al.  is consistent with our model.
With time-resolved imaging, NBR1 localization to individual peroxisomes is expected to be similar to Figure 2B. The signature of cluster formation would be a sudden increase of NBR1 number on peroxisomes with clusters at the same time as NBR1 numbers gradually decrease on peroxisomes without clusters.
(c) Our secondary hypothesis is that p62 inhibits NBR1 cluster formation. If true, we would not expect p62 to have a bimodal distribution since p62 would decorate the freely-diffusing NBR1 (which does not change with cluster formation), but not the NBR1 in clusters (which does). We caution that this lack of bimodality may only apply at the early stages of substrate selection, due to the many cellular roles of p62.
If p62 inhibits NBR1 clustering, then when p62 expression is knocked down by siRNA, there should then be an increase in the number of NBR1 clusters on peroxisomes.
3.6. Important Challenges of Selective Autophagy
We have argued that selective autophagy requires an “all or none” signal to decide whether or not to degrade a substrate. We have proposed a receptor clustering mechanism for this all or none signal for peroxisomes, which leads to a prediction that large substrates are more likely to be selected for degradation. This is consistent with some reports of a preference for degradation of large peroxisomes. However, larger peroxisomes might simply be older and/or more damaged. Distinguishing damage-induced selective autophagy from, e.g., disuse-induced selective autophagy is an interesting challenge. Our suggestion that damage-induced selective autophagy should be less size-selective may help in this regard.
Although our cluster selectivity hypothesis does not uniquely explain size-selectivity, our model does present a working hypothesis for the basic mechanism for substrate selection in selective autophagy. Furthermore, the long lifetime of receptor clusters (see Figure 1C or Figure 3B) allows ample time for downstream processes to positively recognize the all or none signal provided by the receptor clusters, and to avoid “false-positive” triggers on e.g., stochastic fluctuations. Such a stable all-or-none signal is a challenge both for time-resolved microscopy studies of receptor dynamics, but also a challenge for any competing models of the selectivity mechanism of selective autophagy that arise in the future.
4. Computational Methods
4.1. NBR1 and p62 Structure
Each NBR1 molecule contains several regions that are essential for pexophagy. The LC3-interacting region (LIR) interacts with the proteins of the autophagy system [25, 47]. The UBA region can bind to ubiquitin [25, 27], and allows attachment to ubiquitin-tagged substrates . The Phox and Bem1p (PB1) region can bind PB1 regions on other proteins , and coiled-coil regions promote self-interaction . The distinctive “J” region allows NBR1 to anchor to membranes .
Similarly, p62 also has LIR, UBA, and PB1 regions [25, 26, 41, 47]. Distinctively, the PB1 region of p62 can bind two other PB1 regions, forming chains of p62 , unlike NBR1 which can only bind one other PB1 . Unlike NBR1, p62 has no J region.
4.2. NBR1 and p62 Association with Peroxisomes and Each Other
NBR1 and p62 both contain ubiquitin-interacting UBA regions [27, 66]. However, these UBA regions have relatively weak affinities compared to expected cellular abundances. For NBR1, Kd,UBA = 3 − 4μM  while typical abundances in human cell lines are no more than ≈125 ppm , or a concentration [NBR1] ≲ 0.6 μM since 1 ppm corresponds to approximately 5 nM . For p62, Kd,UBA = 540 − 750 μM [38, 39] with typical abundances in human cell lines no more than ≈300 ppm, or [p62] ≲1.5 μM. While phosphorylation of p62 significantly increases polyubiquitin association, the enhancement appears to be no more than three-fold . Phosphorylation decreases the association of NBR1 with ubiquitin . The relatively weak affinities of NBR1 and p62 UBA regions with ubiquitin implies that there should only be a small fraction of these receptors on surface-displayed ubiquitin.
How can NBR1 significantly associate with peroxisomes, if not by association with ubiquitin? The essential J region of NBR1 mediates membrane association even without ubiquitin . Kd of typical amphipathic helices can be as low as 20 nM . Coincidence of ubiquitin and membrane association [19, 72] could further decrease the effective Kd of membrane association through an enhanced on rate. Once freely associated with the peroxisomal membrane, NBR1 molecules can interact through coiled-coil domains [19, 28]. While no specific NBR1-NBR1 interactions have been identified, membrane-bound proteins can form clusters through non-specific interactions. The activity of the holin protein, leading to precise lysis timing  and showing cooperative effects across distinct holin species , is thought to follow from clustering due to non-specific interactions. More generally, non-specific clustering of membrane-bound proteins can result from attractive lipid-mediated protein-protein interactions [30–32]. NBR1 oligomerization is also suggested by experiments that show an important role for NBR1 in the formation of protein aggregates prior to their degradation by autophagy . Given a weak attractive interaction, the physical theory of phase separation  indicates that a sufficiently high concentration of NBR1 on a membrane will lead to a concentrated NBR1 phase—i.e., cluster formation. Such formation of homo-oligomeric clusters following NBR1 membrane association through the J region is our cluster-driven selectivity hypothesis.
While p62 has no identified membrane binding domain, its PB1 region has a strong affinity (Kd = 4 − 10 nM 40) which can lead to association with corresponding NBR1 PB1 regions and also to self-association into p62 filaments [41, 42]. While p62 has two binding faces on its PB1 region, NBR1 only has one  and so cannot form filaments.
We know that knockdown of p62 significantly affects pexophagy , but the mechanism is unknown. Given the weak affinity of p62 to ubiquitin, it appears unlikely that p62 competition for NBR1 binding to ubiquitin is significant. While p62 binding to NBR1 and subsequent polymerization of p62 through PB1 domains would increase the number of LIR domains associated with an organelle, this in itself would be in proportion to the amount of NBR1 associated with the organelle and would not affect NBR1 clustering. However, polymer physics leads us to hypothesize that polymeric chains of p62 associated with NBR1 reduce NBR1 self-association through steric repulsion. Such a repulsion arises as the entropic contribution to the free-energy of polymers is decreased when brought close together . Such steric repulsion of membrane associated proteins can inhibit lipid phase separation , or prevent growth of protein clusters .
Accordingly, we computationally model the association and disassociation of NBR1 on the surfaces of multiple peroxisomes—as recruited by ubiquitin. NBR1 can form clusters when its surface concentration is sufficiently high, and these clusters subsequently grow or shrink as determined by the available NBR1. In the model, p62 can be recruited to membrane-associated NBR1, and subsequently polymerize. Because of steric repulsion, NBR1 that is associated with p62 does not participate in cluster formation or growth.
4.3. Rates of NBR1 and p62 Association
While some equilibrium association constants are determined for NBR1 and p62, kinetic rate constants have not yet been measured. We use diffusion-limited association rates [74, 75]. In our model, we require rates for NBR1 to bind to ubiquitin targets on the surface of a peroxisome, or for p62 to bind to NBR1 targets on the surface of a peroxisome. The diffusion-limited arrival rate is known for arrival at small circular targets on a larger sphere , and within our model it is used for the arrival rate of NBR1 and p62. For NBR1, we use an association rate of NBR1 to each peroxisome
with DNBR1 is the bulk NBR1 diffusivity, ρNBR1 is the bulk NBR1 concentration, R the peroxisome radius, Nub the number of ubiquitin on the peroxisome, and sub the target ubiquitin radius. As discussed in the model motivation, we have sufficiently small bulk concentrations of NBR1 and p62 that a negligible fraction of ubiquitin are occupied, so the number of ubiquitin available for binding remains constant as NBR1 bind. We assume that NBR1 transiently bound to ubiquitin immediately associate with the peroxisomal membrane using their J regions. This allows significant NBR1 to accumulate on the peroxisome.
For p62, we similarly have diffusion-limited rates to surface associated NBR1
where Dp62 is the bulk p62 diffusivity, ρp62 is the bulk p62 concentration, sNBR1 is the (target) NBR1 radius, Nvapour is the total number of NBR1 on the peroxisome surface that are not in clusters, and NNBR1(ℓ) is the number of NBR1 with a p62-chain of length ℓ. Association extends the p62 chain length to ℓ+1.
4.4. NBR1 and p62 Dissociation
NBR1 on the surface of the peroxisome can dissociate from the membrane and return to the cytosol. We also model this as a diffusion-limited process, and for circular targets on a larger sphere, the dissociation rates have been determined . From Ghosh et al. , the effective dissociation rate for NBR1 is
where koff,NBR1 is the dissociation rate of NBR1 from the membrane and γNBR1 ≡ Nubsub/(Nubsub+πR) is the fraction of NBR1 that immediately rebind . By equating the on and off rates of NBR1 to the peroxisome surface, Jon,NBR1 and Joff,NBR1, we see that the steady-state Nvapour is proportional to Nub, and otherwise independent of the peroxisomal radius R:
Equivalently, the steady-state surface concentration of NBR1 is proportional to that of ubiquitin.
In our model we assume that when NBR1 dissociates, any associated p62 chains dissociate as well. In addition, PB1 bonds (p62-p62 or p62-NBR1) within membrane-associated polymers will each break at a rate
where γp62 ≡ NvapoursNBR1/(NvapoursNBR1 + πR) . When a PB1 bond breaks, the portion of the p62 chain beyond the bond (i.e., further from the NBR1 than the bond) dissociates. Since Equation (6) is the rate per PB1 bond, the rate of any PB1 in a chain of length ℓ breaking is proportional to ℓ.
4.5. NBR1 Cluster Formation
Large domains or clusters are generally only thermodynamically stable above some saturation density σ∞ of particles. Smaller clusters are less stable and require a higher density of particles to avoid shrinkage and evaporation; this is known as the Gibbs-Thomson effect . The Gibbs-Thomson effect is due to the increased curvature of the edge of a small cluster compared to a large cluster. Qualitatively, the increased curvature both reduces local bonding and allows particles more directions to escape. In equilibrium, these lead to a higher vapor concentration near the cluster edge to balance the increased escape rate. We need to consider the Gibbs-Thompson effect for small NBR1 clusters.
Before cluster nucleation there will be freely diffusing NBR1 on a peroxisome of area 4πR2 with surface concentration σ and total NBR1 of N = 4πR2σ, where R is the peroxisomal radius. After cluster nucleation, there will be a cluster of radius r with Nclust NBR1 in equilibrium with a surface concentration σgt. Since the cluster is small, the surface concentration will satisfy . Since the number of NBR1 doesn't change, we have N = Nclust + Nsurf. We also satisfy the Gibbs-Thomson effect , with the cluster radius r,
where ν is a constant “capillary length” associated with NBR1 clusters.
We cannot simply allow nucleation when σ ≥ σgt, since the original surface concentration must also provide the NBR1 for the cluster formation. If b is the area per NBR1 in a cluster, then the cluster area . Since , then we require
While we can therefore accommodate a range of possible cluster sizes r, there will be a smallest r determined by minimizing Equation (8) with respect to r. This determines a critical (minimal) surface concentration that allows for cluster nucleation, σ*, where
We see that the minimum supersaturation required for nucleation is lower for larger peroxisomes . We also determine the cluster size after nucleation, .
Following our assumption that NBR1 associated with p62 does not participate in cluster formation, only NBR1 with no associated p62 chain (NNBR1(l = 0)) contribute to the concentration required for nucleation in Equation (9), so that we must have
We have assumed that each peroxisome will harbor either one or zero clusters. For other small biological systems with clusters, including bacterial holin domains [8, 9] and yeast polarity clusters , multiple clusters rapidly resolve to a single cluster. We also note that any supersaturation will be quickly absorbed by the first cluster to nucleate, suppressing further cluster nucleation by Equation (9).
4.6. NBR1 Cluster Growth
Existing NBR1 clusters can gain NBR1 and grow, or lose NBR1 and shrink. To determine the growth of an existing NBR1 cluster on a peroxisome we adapt the derivation in Appendix B of Brown and Rutenberg .
Only NBR1 without a p62 chain (with ℓ = 0) can contribute to cluster growth. The number of such NBR1 on the peroxisome surface, N(0), divided by the surface area 4πR2, gives a surface concentration f0. In principle, f0 is a spatial field over the peroxisome surface, varying depending on location. The dynamics of f0 is then described by the partial differential equation
The first term on the right hand side captures diffusion of NBR1 on the surface, with Ds the surface diffusivity and ∇2 a two-dimensional Laplacian. The second term describes the arrival of NBR1 from the cytosol. The third term represents the dissociation of NBR1 from the surface. The fourth term gives the addition of p62 to an NBR1 with no p62, so that it has p62 and can no longer participate in cluster growth. The fifth term reflects the dissociation of p62 chains from NBR1, which allows those NBR1 to then participate in cluster formation.
After substituting the full expressions for Jon,NBR1 from Equation (2), Joff,NBR1 from Equation (3), Joff,NBR1 from Equation (4), and Joff,p62 from Equation 6 into Equation 11, and assuming close to steady state (i.e., df0/dt = 0) we obtain
where we define ,
Equation (12) is solved in Appendix B of Brown and Rutenberg  to determine the net flux to the cluster. This determines the dynamics of a cluster with Nclust molecules on a peroxisome of radius R,
We use Equation (15) to determine the change in time of the size of every NBR1 cluster in our model.
4.7. Kinetic Model
We implement our kinetic rates to continually update the NBR1 in our system. Approximately 50% of NBR1 colocalizes with catalase and PMP70 , a peroxisome matrix and membrane protein, respectively, indicating that NBR1 dynamics on peroxisomes are probably not significantly buffered by other cellular processes. For every peroxisome, we track each NNBR1(ℓ) and Nclust. The peroxisomal NBR1 is then , and we can sum that over all peroxisomes to obtain Ntot,peroxisomal. We conserve the total amount of NBR1 in the system, so that Ntot,bulk = Ntot − Ntot,peroxisomal and the bulk density ρNBR1 = Ntot,bulk/V, where V is the total system volume.
Only approximately 10% of p62 colocalizes with peroxisomes , which is consistent with the many roles of p62 for autophagy  as well as other cellular pathways . As a result, we expect that uptake by peroxisomes of p62 will not significantly change cytosolic concentrations. Accordingly, we hold p62 concentrations constant.
4.8. Selectivity Calculation
In Figure 3B, we show how the remaining peroxisomal volume fraction depends on size-selectivity when a fixed number fraction of peroxisomes are degraded. With zero selectivity, with targets randomly chosen, the remaining volume and number fractions are proportional. With maximum selectivity, i.e., with only the larger peroxisomes selected for degradation, the volume fraction remaining at number fraction r is
where Rmin is the radius of the smallest peroxisome, Rmax(r) is the radius of the largest peroxisome not selected for degradation, P(R) is the distribution function of peroxisomal radii, and Vi is the total initial peroxisomal volume. The corresponding number fraction is
and we see that Rmax is determined by r.
4.9. Parameters, Initial Conditions, and Numerical Details
The radius of a globular protein is approximately R = 0.066M1/3 for R in nm with M in Daltons . We use this radius to estimate the size of diffusive targets on peroxisomes. For ubiquitin, of mass 8 kDa [79, 80], rub = 1.32 nm. For NBR1, with mass of approximately 107 kDa [29, 81], rNBR1 = 3.14 nm.
While the diffusivity of NBR1 or p62 have not been measured, we can scale the diffusivity of EYFP which is approximately DYFP = 1 μm2/s  (with mass MYFP = 27 kDa). Assuming spherical (globular) proteins and corresponding Stokes-Einstein diffusivity, the diffusivity scales with inverse radius (or cube root of the mass), and we obtain DNBR1 = 0.63 μm2/s and with a p62 mass of 62 kDa [83, 84] Dp62 = 0.83 μm2/s.
Within our model, systems of many peroxisomes have peroxisome radii distributed exponentially, qualitatively like measured peroxisome size distributions [85, 86]. In ensemble systems, , where P(R) is the probability of a peroxisome of radius R, and we use Rs = 0.1 μm. The number of ubiquitin on a given peroxisome will be proportional to the surface area, , with the ubiquitin coefficient n0 typically 100 unless otherwise stated, and R0 = 0.25 μm.
We use a system volume V = Npv, where NP is the number of peroxisomes, and v = 10 μm3 is the volume per peroxisome, unless otherwise stated. 300 peroxisomes has been reported as an average number for mammalian cells . Therefore, the volume inside a spherical cell of radius 10 μm, divided among 300 peroxisomes, is approximately 10 μm3 per peroxisome.
For cluster formation, we assume the capillary length ν is the size of a single NBR1 protein, so ν = rNBR1 = 3.14 nm, and that the area per molecule is b = ν2 = 9.86 nm2. This is consistent with capillary lengths of one  and several  particle widths for 2d and 3d systems, respectively. The vapor pressure σ∞ is taken to be 10 μm−2 on the peroxisome membrane. For a typical peroxisome of radius R = 0.25 μm, this is approximately a single molecule on the peroxisomal surface.
We use koff,NBR1 = 0.1 μm s−1 for the dissociation rate of NBR1 from the peroxisome membrane, which yields a Kd value similar to a 20 nM value measured for amphipathic helices . We choose koff,p62 = 0.04 s−1 for the dissociation rate of p62 from NBR1 and p62 (interaction through the PB1 domain [40, 41, 65]) using the Kd range of 4 - 10 nM for PB1-PB1 bonds .
4.10. Model Limitations
In this section we address some of the limitations of our approach. Our results should not be qualitatively affected by these limitations, whereas precise quantitative predictions would need a more realistic and dynamic cellular geometry, precise parameterization, and a fully stochastic multi-scale approach.
We have used deterministic dynamics in our modeling approach; this was necessary for computational efficiency since our system spans many length and time scales. Our model does not include stochastic effects for the change in molecule number on peroxisomes or in clusters, or for the nucleation of clusters. Change in molecule numbers on peroxisomes and in clusters is determined by the net flux, with frequent molecular association and dissociation. Our deterministic approach reflects an average behavior, and has two limitations. First it does not account for the discrete nature of molecules on peroxisomes or in clusters. The discreteness will affect our nucleation threshold, and we could impose that the number of molecules in a nucleated cluster, , is at least one. Since is strongly R dependent, a larger minimum cluster size due to discreteness will raise the minimum peroxisome size R* that nucleates clusters. This will enhance our predicted size selectivity effect, and so amounts to a conservative approximation. Second, we assume that cluster nucleation occurs deterministically at the threshold. Since nucleation rates typically strongly increase with concentration, the threshold approximation is expected to reproduce the qualitative nucleation behavior.
We have also abstracted the cellular context into uniform concentrations of bulk solutes, rather than an exact stochastic particle-based approach that includes cellular synthesis and degradation—again for computational efficiency. A significant additional advantage of a uniform solute approximation is that it does not require us to model the precise cellular geometry, such as peroxisome locations. This does mean, however, that we cannot treat screening effects between peroxisomes. Identical peroxisomes will behave identically in our model but not within the cell.
An important stochastic effect that we do not include is the downstream autophagy process that removes peroxisomes from the system. We would expect this (missing) process to limit the growth of NBR1 clusters. Without it, NBR1 clusters, according to Equation (15), could grow indefinitely. For the parameters of our model, clusters on typical peroxisomes of radius R = 0.25 μm could reach 10–20% surface coverage at late times. The absence of downstream degradation is most significant for the maximum selectivity case of Figure 3, since in that case the largest remaining peroxisome should retain its cluster until degradation.
We have simplified our receptor dynamics as much as possible. For example, no p62-associated NBR1 will form clusters, but no NBR1 in clusters will associate with p62. We have also assumed that NBR1 is recruited to peroxisomes only by membrane-associated ubiquitin. It is also possible that NBR1 could also be recruited by already membrane-associated NBR1. Such direct recruitment could affect the cluster size dynamics .
As noted previously, many of our parameters are generic estimations, as they have not been measured directly. While the qualitative physics of nucleation and diffusion will be unchanged by large parameter changes, the degree of size-selectivity and the timing and extent of cluster nucleation is parameter dependent. We expect that by not including stochastic effects, we will have effectively shifted the appropriate parameter values. Since our parameters are rough estimations in any case, this does not change our qualitative results and conclusions.
AB implemented the computational model and analyzed the data; AB and AR developed the research direction, developed the computational model, and wrote the paper.
AR thanks the Natural Science and Engineering Research Council of Canada (NSERC) for operating grant support (RGPIN-2014-06245). We thank ACENET for computational resources. AB also thanks NSERC, the Killam Trusts, and the Walter C. Sumner Memorial Foundation for fellowship support.
Conflict of Interest Statement
We thank Peter Kim for helpful discussions.
2. Rogov V, Dotsch V, Johansen T, Kirkin V. Interactions between autophagy receptors and ubiquitin-like proteins form the molecular basis for selective autophagy. Mol Cell. (2014) 53:167–78. doi: 10.1016/j.molcel.2013.12.014
10. Cemma M, Kim PK, Brumell JH. The ubiquitin-binding adaptor proteins p62/SQSTM1 and NDP52 are recruited independently to bacteria-associated microdomains to target Salmonella to the autophagy pathway. Autophagy (2011) 7:341–5. doi: 10.4161/auto.7.3.14046
11. Wild P, Farhan H, McEwan DG, Wagner S, Rogov VV, Brady NR, et al. Phosphorylation of the autophagy receptor optineurin restricts Salmonella growth. Science (2011) 333:228–33. doi: 10.1126/science.1205405
12. Mostowy S, Sancho-Shimizu V, Hamon MA, Simeone R, Brosch R, Johansen T, et al. p62 and NDP52 proteins target intracytosolic Shigella and Listeria to different autophagy pathways. J Biol Chem. (2011) 286:26987–95. doi: 10.1074/jbc.M111.223610
17. Monastryska I, Sjollema K, van der Klei IJ, Kiel JAKW, Veenhis M. Microautophagy and macropexophagy may occur simultaneously in Hansenula polymorpha. FEBS Lett. (2004) 568:135–8. doi: 10.1016/j.febslet.2004.05.018
18. Kim PK, Hailey DW, Mullen RT, Lippincott-Schwartz J. Ubiquitin signals autophagic degradation of cytosolic proteins and peroxisomes. Proc Natl Acad Sci USA. (2008) 105:20567–74. doi: 10.1073/pnas.0810611105
22. Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, Sethna JP. Universally sloppy parameter sensitivities in systems biology models. PLoS Comput Biol. (2007) 3:e30189. doi: 10.1371/journal.pcbi.0030189
27. Vadlamudi RK, Joung I, Strominger JL, Shin J. p62, a phosphotyrosine-independent ligand of the SH2 domain of p56, belongs to a new class of ubiquitin-binding proteins. J Biol Chem. (1996) 271:20235–7. doi: 10.1074/jbc.271.34.20235
28. Kirkin V, Lamark T, Sou YS, Bjorkoy G, Nunn JL, Bruun JA, et al. A role for NBR1 in autophagosomal degradation of ubiquitinated substrates. Mol Cell. (2009) 33:505–16. doi: 10.1016/j.molcel.2009.01.020
29. Nicot AS, Verso FL, Ratti F, Pilot-Storck F, Streichenberger N, Sandri M, et al. Phosphorylation of NBR1 by GSK3 modulates protein aggregation. Autophagy (2014) 10:1036–53. doi: 10.4161/auto.28479
30. Heimburg T, Biltonen RL. A Monte Carlo simulation study of protein-induced heat capacity changed and lipid-induced protein clustering. Biophys J. (1996) 70:84–96. doi: 10.1016/S0006-3495(96)79551-6
31. Gil T, Ipsen JH, Mouritsen OG, Sabra MC, Sperotto MM, Zuckermann MJ. Theoretical analysis of protein organization in lipid membranes. Biochim Biophys Acta. (1998) 1376:245–66. doi: 10.1016/S0304-4157(98)00022-7
32. Lague P, Zuckermann MJ, Roux B. Lipid-mediated interactions between intrinsic membrane proteins: dependence on protein size and lipid composition. Biophys J. (2001) 81:276–84. doi: 10.1016/S0006-3495(01)75698-6
38. Long J, Gallagher TRA, Cavey JR, Sheppard PW, Ralston SH, Layfield R, et al. Ubiquitin recognition by the ubiquitin-associated domain of p62 involves a novel conformational switch. J Biol Chem. (2008) 283:5427–40. doi: 10.1074/jbc.M704973200
40. Wilson MI, Gill DJ, Perisic O, Quinn MT, Williams RL. PB1 domain-mediated heterodimerization in NADPH oxidase and signaling complexes of atypical protein kinase C with Par6 and p62. Mol Cell. (2003) 12:39–50. doi: 10.1016/S1097-2765(03)00246-6
41. Lamark T, Perander M, Outzen H, Kristiansen K, Overvatn A, Michaelsen E, et al. Interaction codes within the family of mammalian Phox and Bem1p domain-containing proteins. J Biol Chem. (2003) 278:34568–81. doi: 10.1074/jbc.M303221200
43. Hristova K, Needham D. The influence of polymer-grafted lipids on the physical properties of lipid bilayers: a theoretical study. J Colloid Inter Sci. (1994) 168:302–14. doi: 10.1006/jcis.1994.1424
44. Sieber JJ, Willig KI, Kutzner C, Gerding-Reimers C, Harke B, Donnert G, et al. Anatomy and dynamics of a supramolecular membrane protein cluster. Science (2007) 317:1072–6. doi: 10.1126/science.1141727
48. Veenhuis M, Zwart K, Harder W. Degradation of peroxisomes after transfer of methanol-grown Hansenula polymorpha into glucose-containing media. FEMS Microbiol Lett. (1978) 3:21–28. doi: 10.1111/j.1574-6968.1978.tb01876.x
52. Veenhuis M, Douma A, Harder W, Osumi M. Degradation and turnover of peroxisomes in the yeast Hansenula polymorpha induced by selective inactivation of peroxisomal enzymes. Arch Microbiol. (1983) 134:193–203. doi: 10.1007/BF00407757
53. Vagne Q, Turner MS, Sens P. Sensing size through clustering in non-equilibrium membranes and the control of membrane-bound enzymatic reactions. PLoS ONE (2015) 10:e0143470. doi: 10.1371/journal.pone.0143470
54. Odendall C, Dixit E, Stavru F, Bierne H, Franz KM, Durbin AF, et al. Diverse intracellular pathogens activate type III interferon expression from peroxisomes. Nat Immunol. (2014) 15:717–26. doi: 10.1038/ni.2915
55. Hou F, Sun L, Zheng H, Skaug B, Jiang QX, Chen ZJ. MAVS forms functional prion-like aggregates to activate and propagate antiviral innate immune response. Cell (2011) 146:448–61. doi: 10.1016/j.cell.2011.06.041
56. Onoguchi K, Onomoto K, Takamatsu S, Jogi M, Takemura A, Morimoto S, et al. Virus-infection or 5′ppp- RNA activates antiviral signal through redistribution of IPS-1 mediated by MFN1. PLoS Pathog. (2010) 6:e1001012. doi: 10.1371/journal.ppat.1001012
59. Kiel JAKW, Emmrich K, Meyer H, Kunau W. Ubiquitination of the the peroxisomal targeting signal type 1 receptor, Pex5p, suggests the presence of a quality control mechanism during peroxisomal matrix protein import. J Biol Chem. (2005) 280:1921–930. doi: 10.1074/jbc.M403632200
62. Shibata M, Oikawa K, Yoshimoto K, Kondo M, Mano S, Yamada K, et al. Highly oxidized peroxisomes are selectively degraded via autophagy in Arabidopsis. Plant Cell. (2013) 25:4967–83. doi: 10.1105/tpc.113.116947
63. van Zutphen T, Veenhuis M, van der Klei IJ. Damaged peroxisomes are subject to rapid autophagic degradation in the yeast Hansenul polymorpha. Autophagy (2011) 7:863–72. doi: 10.4161/auto.7.8.15697
65. Ciuffa R, Lamark T, Tarafder AK, Guesdon A, Rybina S, Hagen WJH, et al. The selective autophagy receptor p62 forms a flexible filamentous helical scaffold. Cell Rep. (2015) 11:748–758. doi: 10.1016/j.celrep.2015.03.062
66. Walinda E, Morimoto D, Sugase K, Konuma T, Tochio H, Shirakawa M. Solution structure of the ubiquitin-associated (UBA) domain of human autophagy receptor NBR1 and its interaction with ubiquitin and polyubiquitin. J Biol Chem. (2014) 289:13890–902. doi: 10.1074/jbc.M114.555441
67. Geiger T, Wehner A, Schaab C, Cox J, Mann M. Comparative proteomic analysis of eleven common cell lines reveals ubiquitous but varying expression of most proteins. Mol Cell Proteomics. (2012) 11:M111.014050. doi: 10.1074/mcp.M111.014050
69. Matsumoto G, Wada K, Okuno M, Kurosawa M, Nukina N. Serine 403 phosphorylation of p62/SQSTM1 regulates selective autophagic clearance of ubiquitinated proteins. Mol Cell. (2011) 44:279–89. doi: 10.1016/j.molcel.2011.07.039
70. Mardakheh FK, Auciello G, Dafforn TR, Rappoport JZ, Heath JK. Nbr1 is a novel inhibitor of ligand-mediated receptor tyrosine kinase degradation. Mol Cell Biol. (2010) 30:5672–85. doi: 10.1128/MCB.00878-10
71. Stahelin RV, Long F, Peter BJ, Murray D, Camilli PD, McMahon HT, et al. Contrasting membrane interaction mechanisms of AP180 N-terminal homology (ANTH) and epsin N-terminal homology (ENTH) domains. J Biol Chem. (2003) 278:28993–9. doi: 10.1074/jbc.M302865200
78. Erickson HP. Size and shape of protein molecules at the nanometer level determined by sedimentation, gel filtration, and electron microscopy. Biol Proced Online (2009) 11:32–51. doi: 10.1007/s12575-009-9008-x
80. van Delft S, Govers R, Strous GJ, Verkleij AJ, van Bergen PMP, En Henegouwen B. Epidermal growth factor induces ubiquitination of Eps15. J Biol Chem. (1997) 272:14013–6. doi: 10.1074/jbc.272.22.14013
81. Waters S, Marchbank K, Solomon E, Whitehouse C, Gautel M. Interactions with LC3 and polyubiquitin chains link nbr1 to autophagic protein turnover. FEBS Lett. (2009) 583:1846–52. doi: 10.1016/j.febslet.2009.04.049
83. Pankiv S, Clausen TH, Lamark T, Brech A, Bruun JA, Outzen H, et al. p62/SQSTM1 binds directly to Atg8/LC3 to facilitate degradation of ubiquitinated protein aggregates by autophagy. J Biol Chem. (2007) 282:24131–45. doi: 10.1074/jbc.M702824200
85. Liu F, Lu Y, Pieuchot L, Dhavale T, Jedd G. Import oligomers induce positive feedback to promote peroxisome differentiation and control organelle abundance. Dev Cell. (2011) 21:457–68. doi: 10.1016/j.devcel.2011.08.004
86. Vizeacoumar FJ, Torres-Guzman JC, Bouard D, Aitchison JD, Rachubinski RA. Pex30p, Pex31p, and Pex32p form a family of peroxisomal integral membrane proteins regulating peroxisome size and number in Saccharomyces cerevisiae. Mol Biol Cell. (2004) 15:665–77. doi: 10.1091/mbc.E03-09-0681
88. Strobel M, Heinig KH, Moller W. Three-dimensional domain growth on the size scale of the capillary length: effective growth exponent and comparative atomistic and mean-field simulations. Phys Rev B. (2001) 64:245422. doi: 10.1103/PhysRevB.64.245422
Keywords: pexophagy, selective autophagy, peroxisomes, NBR1, p62, receptor clustering, biological physics, computational modeling
Citation: Brown AI and Rutenberg AD (2017) A Model of Autophagy Size Selectivity by Receptor Clustering on Peroxisomes. Front. Phys. 5:14. doi: 10.3389/fphy.2017.00014
Received: 14 February 2017; Accepted: 02 May 2017;
Published: 19 May 2017.
Edited by:David Holcman, École Normale Supérieure, France
Reviewed by:Nir Shachna Gov, Weizmann Institute of Science, Israel
Luis Diambra, National University of La Plata, Argentina
Copyright © 2017 Brown and Rutenberg. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.