A game-theoretic approach to deciphering the dynamics of amyloid-β aggregation along competing pathways

Aggregation of amyloid-β (Aβ) peptides is a significant event that underpins Alzheimer's disease (AD). Aβ aggregates, especially the low-molecular weight oligomers, are the primary toxic agents in AD pathogenesis. Therefore, there is increasing interest in understanding their formation and behaviour. In this paper, we use our previously established results on heterotypic interactions between Aβ and fatty acids (FAs) to investigate off-pathway aggregation under the control of FA concentrations to develop a mathematical framework that captures the mechanism. Our framework to define and simulate the competing on- and off-pathways of Aβ aggregation is based on the principles of game theory. Together with detailed simulations and biophysical experiments, our models describe the dynamics involved in the mechanisms of Aβ aggregation in the presence of FAs to adopt multiple pathways. Specifically, our reduced-order computations indicate that the emergence of off- or on-pathway aggregates are tightly controlled by a narrow set of rate constants, and one could alter such parameters to populate a particular oligomeric species. These models agree with the detailed simulations and experimental data on using FA as a heterotypic partner to modulate the temporal parameters. Predicting spatio-temporal landscape along competing pathways for a given heterotypic partner such as lipids is a first step towards simulating scenarios in which the generation of specific ‘conformer strains’ of Aβ could be predicted. This approach could be significant in deciphering the mechanisms of amyloid aggregation and strain generation, which are ubiquitously observed in many neurodegenerative diseases.


Introduction
Aggregation of the protein amyloid-β (Aβ) is one of the central processes in the aetiology of Alzheimer's disease (AD). Generated by the proteolytic processing of amyloid precursor protein (APP), Aβ peptides (Aβ40 or Aβ42) spontaneously aggregate to form insoluble fibrils that deposit as senile plaques in the AD brain. During aggregation, the low-molecular weight oligomers formed are known to be the primary toxic species responsible for synaptic dysfunction and neuronal loss [1][2][3][4][5][6]. An increasing number of reports indicate that structural polymorphism and heterogeneity within the aggregates could contribute to clinical phenotypes observed among AD patients [7,8]. Therefore, over the last decade, significant efforts are focused on understanding the biophysical and biochemical aspects of aggregation as well as the molecular understanding of the aggregates.
Aβ aggregation follows a nucleation-dependent, sigmoidal growth kinetics involving a key ratelimiting event of nucleus or nuclei formation [9][10][11][12][13]. Since the nucleation plays an important role in determining the morphology of the fibrils formed, the dynamics associated with reactions leading up to nucleation are critical determinants of aggregation. In this regard, the sensitivity of Aβ to environmental factors and many interacting partners due to its intrinsic disorder and amphipathic nature [14][15][16][17][18] play a key role in Aβ adopting multiple pathways depending on the aggregation conditions. An important ramification of this is that the oligomers may not be obligate intermediates of fibril formation but those with distinct conformations can be formed along alternative aggregation pathways (off-pathways) [13,[19][20][21][22][23]. This is significant because such interactions, depending on the structure of the oligomer, determine the morphology of the aggregates formed and consequently, the toxicity and phenotypes.
Therefore, it is imperative to gain an understanding of how physiological interacting partners of Aβ affect its aggregation dynamics. Being generated from the membrane-spanning domain of the APP, Aβ displays synchronous and perpetual interaction with membrane lipids [24][25][26][27][28][29][30]. Interfaces of lipids and fatty acids (FAs) are also abundant in both cerebral vasculature and cerebral spinal fluid (CSF) [31,32]. Previous reports have established that phase transitions of surfactants and membrane lipids modulate Aβ aggregation in a concentration-dependent manner to generate aggregates by an alternative, offpathway from the canonical fibril formation, on-pathway [13,16,20,[33][34][35][36][37]. Specifically, in micellar lipids, low-molecular weight oligomers were generated along off-pathway in the presence of fatty acid near and above their respective critical micelle concentrations (CMC) ( pseudo-micellar and micellar, respectively) and not below CMC (non-micellar) which augmented the fibril formation in the onpathway [16,34,38].
The modulation of aggregation by heterotypic interactions between Aβ and lipids posit the question of what spatio-temporal parameters govern the modulatory dynamics, and whether one could simulate the temporal emergence and disappearance of aggregates as a function of heterotypic Aβ-lipid interactions. In this work, we have approached to answer these questions using a competition-based (built qualitatively upon the idea of game theory) approach to determine the dynamics in the temporal evolution of Aβ aggregates along the pathways influenced by fatty acid surfactants (L). Our rationale for such an approach is that the stochasticity and the often exclusive pathways of Aβ aggregation present 'win or loss' scenarios with respect to pathway adoption, tightly governed by the concentration and phase transitions of L. The mathematical analysis of this problem was taken up in two layers, one feeding into the other. The first is a six species, coarse-grained, reduced-order model (ROM), while the second is a more detailed model called ensemble kinetic simulation (EKS), which captures the temporal kinetics of reactions at the atomistic scale (considered as point particles). The ROM approach lends itself to a detailed analysis in a manner that cannot be performed in highresolution models as we have shown before [34,39]. Phenomenologically inspired by the biophysical framework, and using toy models, ROM provides insights into the dynamics of mechanism that are previously unknown. In addition, the outcomes of the ROM analysis provide the appropriate cues to investigate the mechanism deeper with the EKS models. These models are partly validated by bulk kinetic and thermodynamic features using biophysical experiments. The simulations, supported by biophysical analyses, provide a temporal contour map along competing pathways, and present a unique perspective on otherwise unknown evolution of aggregates along multiple pathways.

Reduced-order kinetic modelling
The model presented here consists of a reduced order, comprising only six species of Aβ that interact with the fatty acid surfactant, L. Even with just six species, there are infinitely many rate regimes, most of which would be physically inconsequential. Thus, only physically meaningful rate regimes suggested from experiments and our previous studies [34,39] were chosen, and key parameters were varied to understand the dynamics. Specifically, two models were considered: (i) the base model, where the forward rate constants to back constants were taken to be 1000, and (ii) a second 'pathological' model, where the forward and backward reactions are taken to be identical. The second model has no known physical basis; however, it can be considered as a sort of parameter sensitivity study and an extreme case when the physiological process breaks down.
A schematic of such a model is presented below (see also figure 1). In this model, Aβ monomers react with the pseudo-micellar fatty acid surfactants, L to modulate the formation of on-or off-pathway aggregates. The system of chemical reactions in our model consists of the following: The non-prime species, A 1 , A n and A m represent on-pathway Aβ monomers (A 1 ) and oligomers (A n and A m where m is an integral multiple of n); whereas the prime species, A 0 1 , A 0 n and A 0 m , are the corresponding off-pathway species which are generated through a reaction with the pseudo-micellar surfactant, L. The rate constants k + i (i = 1-6) are indicated in the reaction schematic above where the '+' represents a forward rate and '−', a backward rate. These reactions were formulated based on experimental evidence demonstrated earlier [40]. In the computations to follow, for each species, n = 4 and m = 20 unless otherwise specified, which denotes the order of oligomer [33]. The n, m values in the computations were kept low to minimize computational time. This is also because only significant qualitative features in the system were being sought by ROM, and a more fine-grained approach by Figure 1. Schematic of on-and off-pathway aggregation model based on the six-species reaction scheme described earlier.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191814 EKS modelling provides atomistic temporal analyses using the output from ROM. However, it must be noted that the key results of the study were examined for different values of n and m to ensure qualitative similarities and with no loss of generality as shown previously [34]. The reaction scheme was used to develop the corresponding kinetic model comprising a system of six nonlinear differential equations. This system was then put into non-dimensional form. Using A 0 as the characteristic concentration of monomers and 1=k À 1 the characteristic time, the dimensionless species are defined as follows: The reaction constants are similarly defined as follows: Note that both α 3 and α 4 have a factor L which is responsible for off-pathway aggregation. These two parameters serve as the bridge variables between on-and off-pathway species. Using the law of mass action kinetics, the dimensionless system of differential equations was formulated as follows: As stated earlier, primarily two models referred to as the Base Model and Model 2 were analysed, which are distinguished by the choice of fixed parameter values; i.e. the rate constant ratios in the pure on-and off-pathways. In the Base Model, all forward rates (α 1 , α 2 , α 5 , α 6 ) and all backward rates (β 1 , β 2 , β 3 , β 4 , β 5 ) were set to 1 and 0.001 based on previous mathematical models and experimental data [39,40]. In the context of the Base Model, a forward rate is defined as one that converts a smaller oligomer into a larger aggregate, and backward being the reverse process. It must be noted that since ROM is a bulk averaged model, precise one-to-one mapping of its rate constants to that of the detailed EKS model is neither practical nor meaningful. In Model 2, all forward and backward rates were set to unity. Ode 45 solver (Matlab) was used for our numerical computations.
A convenient approach to the problem would be to analyse the model equations (2.1)-(2.6) from a game-theoretic point of view. Such an approach warrants finding the conditions under which the triplet (B 1 , B n , B m ) are greater, less or equal to (B 0 1 , B 0 n , B 0 m ) respectively; equality would indicate the Nash equilibrium. A similar game-theoretic treatment was applied to a simpler system in our earlier work on multiple-pathway protein aggregation [34], and also by others on various biochemical systems [41][42][43]. In the context of amyloid protein aggregation, the current model system shows the emergence of new states discussed in detail in §3.1.3, which have previously not been observed and lead to new experimental questions about dominant chemical reaction fluxes in competing systems.

Ensemble kinetic simulation
Detailed insights into the switching behaviour between on-and off-pathways were formulated by a combined off-on-pathway EKS model. EKS model has previously been applied for Aβ aggregation system [11,12,34,39,[44][45][46]. In this paper, we have extended our previous work by adding switching royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191814 reactions considering off-to-on and on-to-off oligomer conversion. It has to be borne in mind that the switching reactions only take effect from perturbation events such as changes in the concentrations of pseudo-micellar fatty acid, L [16].
In the EKS model, a set of reactions was considered to represent the on-pathway, off-pathway and their switching, and the flux for each reaction was computed. The system of differential equations of each species present in the reaction system were identified and solved using the ODE 23s solver (Matlab). The following is the reaction scheme considered (corresponding differential equations are presented in Appendix C): I. Reactions of on-pathway: (considering A 12 as F): II. Reactions of off-pathway model: III. On-to-off switching reaction: The corresponding flux for the reactions is given as follows: I. On-pathway reactions flux: II. Off-pathway reactions flux: III. Switching flux: Here, and in Appendix C, A i denotes an on-pathway i-mer, A 0 i denotes an off-pathway i-mer, L denotes pseudo-micellar surfactants, F denotes post-nucleation on-pathway aggregates (here A 12 is considered equivalent to F which corresponds to an on-pathway nucleus of 12mer as previously reported [40]; F for the sake of simplicity), F 0 i is an off-pathway oligomer, signal is the total thioflavin-T (ThT) fluorescence intensity which is expressed as the sum of the on-pathway ThT signal (signal on ) and the off-pathway ThT signal (signal off ) (as shown in Appendix B; this uses an arbitrary mapping constant to map the total oligomer concentration to the experimentally observed ThT signal intensity).
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191814 Note that in the EKS models, we consider the most general case where there can be switching between any on-or off-pathway oligomer of size A 1 to A 11 . Similarly, smaller off-pathway oligomers from A 0 16 -A 0 23 were considered as F 0 1 and larger off-pathway oligomers were considered as F 0 i ði ¼ 1, . . . ,4Þ while a dissociation of F 0 4 was considered to lead to the formation of F 00 1 , which is a kinetically trapped offpathway oligomer that does not aggregate further. The existence of such on-and off-pathway oligomers and the validity of our combined on-and off-pathway model (barring the switching reaction) have already been established in earlier work [34,47].

Biophysical analysis
Synthetic, wild-type Aβ42 procured from both Peptide 2.0 and Dr Chaterjee's laboratory at the University of Mississippi was used in this study. ThT, sodium dodecyl sulfate (SDS) and lauric acid (C 12:0) was purchased from Sigma-Aldrich (St Louis, MO). Monoclonal Ab9 or Ab5 antibodies were obtained from the University of Florida Center for Translational Research in Neurodegenerative Diseases.

Protein preparations
Preparation of Aβ42 monomers: Aβ42 peptide (1-1.5 mg) that was kept desiccated at −20°C was dissolved in 50 mM NaOH and was allowed to stand at room temperature for up to 45 min. The dissolved peptide was then fractionated on a Superdex-75 HR 10/30 size exclusion chromatography (SEC) column (GE Life Sciences) on a BIORAD FPLC system that was pre-equilibrated with 20 mM Tris at pH 8.00, to separate any preformed aggregates as previously reported [41]. Fractions were collected at a flow rate of 0.5 ml min −1 and stored at 4°C and were used within 24 h to avoid reaggregation. The concentration of the monomeric fractions was calculated using a Cary 50 UV-Vis spectrophotometer (Varian Inc.). The molar extinction coefficient of 1450 cm −1 M −1 at 276 nm was used (www.expasy.org).

ThT fluorescence aggregation assay
For on-to-off-pathway switching reactions, to 150 µl, 50 µM Aβ reactions incubated in buffer alone, a 50 µM ThT solution in the same buffer was added and fluorescence emission (λ = 482 nm) was collected using microplate reader (BioTek Synergy Microplate Reader) at 37°C using an excitation at 452 nm. A 5 mM sodium laurate (C12 fatty acid) sample pre-equilibrated with 50 µM ThT was added to the reactions at 3, 8 and 24 h to initiate switching of pathways. The data were collected at 10 min time intervals. For off-to-on-pathway switching reactions, the 150 µl, 50 µM Aβ reactions pre-incubated in the presence of 5 mM sodium laurate were diluted 5-or 10-fold at 5 and 10 h using buffered 50 µM Aβ monomers and 50 µM ThT such that only the fatty acid concentration is dropped below its CMC. Appropriate blank reactions were monitored simultaneously and were corrected before data processing.

SDS-PAGE and immunoblotting
Aliquots of the reactions were mixed with sample buffer comprising 1% SDS (1× Laemmli sample buffer) and loaded on a precast 4-12% Bio-Rad gel. For calibration, pre-stained molecular weight markers (Invitrogen Inc.) were used. The gels were then electroblotted on 0.45 µm nitrocellulose membrane (GE Life Sciences). The blots were then heated in the microwave for 1 min and were blocked with 5% non-fat dry milk solution with 1% Tween-20 in PBS for 1.5 h. Subsequently, the blots were probed with monoclonal antibody Ab5 or Ab9 (1 : 1000-1 : 2500 dilutions) which bind to residue 1-16 of Aβ. Anti-mouse horseradish peroxide was added to the blot and the blot was developed with ECL reagent (Thermo Fisher Scientific) and imaged with a Bio-Rad Gel Doc system.

. Steady states
In order to study the stability of the system, bulk rate constants obtained from experiments were used to determine steady state, or concentration at a given time-point [39]. One is especially interested in the nonzero terminal states of each species. Numerical computations indicate that concentration of species continues to change over time for our models, heading towards the steady state. However, in all cases these changes were within 0.1% of previous levels for t greater than some critical time, which was considered acceptable as equilibrium. The equilibrium values were also confirmed through Matlab's fsolve function. In all ROM computations discussed in this paper, the initial conditions were taken to be B 1 (0) = 1 with all other species set initially at zero. As seen from figure 2a and b, as time increased, the concentration levels exhibited asymptotic behaviour and each species eventually achieved equilibrium. The time to reach this steady state was sensitive to the choice of rate constants; the Base Model took longer to reach steady state than the Model 2. Also, in all cases analysed the fibril concentrations B m and B 0 m took the longest to reach equilibrium. Due to low forward rates, the concentration size of B 1 stayed high and stable throughout, but large percentage changes in the concentration of B 0 m were observed periodically. Analysis of the concentration patterns of both B 1 and B 0 m over this period revealed that their growth and decline patterns are a reciprocal image of one another: periods of increase in B 1 were accompanied by decline in B 0 m . From the equilibrium analysis of these models, it was discovered that when the ratio of backward to forward rates is close to 1, the model settles at equilibrium more quickly than when the ratio is large (figure 2c). A power-law regression indicates that time to equilibrium (indicated by t eq ) varies with the ratio of forward to backward rates (r) according to t eq ∝ r −0.615 .  The impact of varying n and m in both models was also investigated. In these cases, increasing n and m increased t eq but yielded similar qualitative results, some of which are also discussed in Appendices A and B. It appears that the larger the oligomer size, the higher the power of the nonlinear terms in the governing equations, the greater the potential for over-and under-shoot as the model evolves over time. Thus, it takes longer to achieve equilibrium. However, Model 2 does not show an increase in t eq as noted earlier.

Bridge parameters
The key parameters in our model are α 3 and α 4 which are referred to as 'bridge' or control parameters, and they govern the reaction dynamics between on-and off-pathway. The effect of varying both on species formation was verified while holding all other reaction rates constant (figure 3). When increasing both α 3 and α 4 , a direct increase in the ratio of off-pathway species to on-pathway species was observed. Since B 0 m =B m is not directly governed by the bridge variables, it was slower to react to changes along the bridges, but eventually exhibited what appears to be exponential growth at higher values of the bridge variables (figure 3a). This is probably due to the fact that B 0 m formation is dependent upon α 3 and α 4 , so that increasing α 3 and α 4 eventually impacts B 0 m . Figure 3a and b underscores the importance of the bridge variables. Interestingly, if α 4 was left unchanged and α 3 was increased, there was limited flow-through from B 0 1 to B 0 n and B 0 m : their ratios to the non-prime species increased slightly above unity, but ceased to grow from there on even as α 3 continued to increase. Therefore, this suggests that the bridge reaction B n $ B 0 n is critical in the formation of the larger oligomers, i.e. the n and m species. The ROM modelling, therefore, reveals that bridges between larger oligomers are more significant than the ones across monomers in terms of promoting off-pathway fibril formation. Additional tests were performed to verify conditions for any species to outperform others by appropriate choice of the rate constants. Forcing B 1 to outperform, for instance, is just a matter of reducing or shutting off all the forward reactions. For species further down the reaction-network, forward reactions were required to increase to get the desired outperformance. In the case of B 0 m , out-performance of this species was obtained in absolute terms by increasing the forward reaction rates α 3 , β 2 and α 6 by an order 10 4 . Out-performance by B m could also be achieved in a similar manner. Such an exercise can be significant in helping to identify pathogenic aggregates and shows the robustness of the network under standard reaction rates. Figure 4d provides a schematic of the four equilibrium pathways that our model can achieve, each sensitive to the choices of parameters α 3 and α 4 . In figure 4a, the first schematic highlighted in red is strictly onpathway, where the non-prime species 'win'. The next highlighted in blue is strictly off-pathway, where all off-pathway species wins. The paths indicated in yellow and green are a mixture of on/off-pathway. Figure 4d depicts the network graph corresponding to each of the phases. Computations were conducted  For the Base Model, (α 3 , α 4 ) = (1, 1); is a critical point at which the concentrations of on-and offpathway species are equal. As α 3 and α 4 were varied, dominance of one set of species or pathways over another emerged. Notable too is the fact that the boundaries between the different equilibrium states were almost linear: the line α 4 = 1 determines the switching between on-and off-pathway dominance of n and m species, and the line α 3 = 1 determines the switching between on-and offpathway dominance of monomers. The table in figure 4c shows the equilibrium states as a function of α 3 and α 4 in the form of a pay-off-like matrix. The Nash equilibrium lies at the point where on-pathway species concentrations are equal to off-pathway species concentrations.

A 'game-theoretic' approach to understanding pathways
A similar computation was performed for Model 2 (figure 4b).
Here too, (α 3 , α 4 ) = (1, 1) was a critical point; however, unlike in the Base Model, it does not strictly define out-performance of B 1 over B 0 1 and vice versa; still B 1 outperforming B 0 1 was seen for α 3 > 1 and low α 4 , and B 0 1 outperforming B 1 for α 3 < 1 and high α 4 . The major difference is that the red and blue regions representing the only on-pathway and only off-pathway, respectively, increased at the expense of green and yellow (mixed pathways).
More importantly, the range of points over which on-pathway wins got bigger when backward rates were on par with forward rates. For α 3 < 1, by increasing β 1 and β 5 , B 1 ! B 0 1 bridge towards off-pathway aggregation was effectively shut off, hence the increase in red in the upper left of figure 4b. The reverse happened in the lower right for α 4 < 1 as we observed greater off-pathway aggregation up to a point. Despite increasing α 3 above 2, upon reducing α 4 , a thin band of dominance of B n and B m over their respective primed off-pathway species was continued to be observed. Once again, this shows that the B n ! B 0 n bridge is more critical for off-pathway aggregation of n and m species than the monomerbridge. If α 4 was reduced, an on-pathway dominance of n and m species even for high α 3 was still obtained. Thus, it is difficult to control the off-pathway aggregation of n and m species by tweaking the B 1 ! B 0 1 bridge.   royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191814

Biophysical evidence for the switching of aggregation pathways
The effect of FAs on Aβ42 aggregation has been well established in the Rangachari laboratory [13,44,45]. Specifically, using sodium laurate at concentrations near and well above its CMC, the surfactant was able to modulate Aβ42 aggregation toward off-fibril formation pathway that was populated by low-molecular weight oligomers. At concentrations well below CMC, the fatty acid adopted an on-fibril formation pathway [16]. To experimentally assess the switching of pathways from on-to off-pathway and vice versa by modulating L concentrations, kinetic rate differences in aggregation was investigated using ThT dye. Switching of on-to off-pathway (depicted schematically in figure 5a) was initiated by the addition of 5 mM C12 FA to 25 µM Aβ42 buffered in 20 mM Tris, 50 mM NaCl at pH 8.0. The addition of C12 FA resulted in an increase in ThT fluorescence without any observable lag time (black square; figure 5b). By contrast, Aβ42 in the absence of C12 FA showed a lag phase of approximately 50 h before an increase in ThT fluorescence was observed (black diamond; figure 5b). This behaviour in the presence of C12 FA has been previously observed to generate 12-24mer oligomers of Aβ along the off-fibril formation pathway [48]. In order to evaluate the propensity of bridging from on-to off-pathway, 5 mM C12 FA was added to the control Aβ42 reaction after 0 h ( positive control (black square), 3 h (red circle), 8 h (blue triangle) and 24 h ( purple inverted triangle). Each of such incubations resulted in an exponential increase in ThT fluorescence suggesting switching of pathways from on to off (figure 5b). Analysis of these samples was also performed using a partially denaturing gel electrophoresis (low SDS and no boiling) and immunoblotting (figure 5b). Injections of C12 FA at 3 and 8 h show the presence of 48-60 kDa band corresponding to 12mer oligomers (lanes 1 and 2, respectively) as compared with the corresponding controls generated upon adding buffer in place of royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191814 C12 FA (lanes 3 and 4), which show monomers and some on-pathway aggregates. This suggests that offpathway oligomers are generated (figure 5c). Similarly, FA injected after 24 h and its corresponding control are shown in lanes 5 and 6, respectively, which shows even after 24 h, C12 FA is able to induce the formation of oligomers to a certain extent, with clearly observable emergence of some onpathway fibrils. These results parallel those observed by ThT fluorescence (figure 5b).
To further quantify the extent of bridging, the aggregates generated after the 24 h injection of C12 FA (or buffer for the control) were fractionated by SEC, after an additional 24 h incubation (figure 5d). Prior to fractionation, the samples were centrifuged at 18 000g for 20 min to remove any high molecular weight fibrils, and the supernatant was loaded on to the column. After 24 h, the control in the absence of C12 FA shows a small peak near the void volume at fraction 17 and a monomer peak at fraction 24 (blue; figure 5d). Fractionation of the control reaction at 48 h (after injection of buffer at 24 h) showed a diminished peak at fraction 17 and a reduced monomer peak at fraction 24 (red; figure 5d). The reduction in the monomer peak correlates to being consumed during aggregation. A similar reduction in the aggregate peak between 24 and 48 h can be explained by the fact that many have formed fibrils that are centrifuged out. On the other hand, fractionation of the sample after 48 h with the injection of 5 mM C12 FA at 24 h, showed a larger peak at fraction 18 and a reduced monomer peak at fraction 25 (black; figure 5d ). This suggests two possibilities: (i) the unreacted monomers adopt off-pathway upon introduction of C12 FA, and/or (ii) the pre-formed aggregates along the on-pathway are rerouted back through the off-pathway, in other words, switching. More detailed analysis on this is discussed later in the article.
To assess a similar switching of pathways, we performed the off-to on-pathway (schematically depicted in figure 5e) switching again using the established C12 FA kinetics. Incubation of 5 mM C12 FA shows an exponential increase in ThT fluorescence (black diamond; figure 5f ). To effect switching of off-to on-pathway after certain time periods, the sample was diluted 5-and 10-fold such that the effective concentration of C12 FA drops to 1 and 0.1 mM, which are well below the CMC of the surfactant. It is well established that well below CMC, Aβ aggregation is augmented [16], and therefore, dilutions of 5 mM C12 FA must induce faster rates of aggregation. When dilutions were introduced, at 5 and 24 h time points (arrows; figure 5f ), appropriately blank subtracted data showed an increase in ThT fluorescence as expected for both dilutions suggesting the switching of off-to on-pathway (figure 5f ). Partially denaturing gel electrophoresis and immunoblotting further confirmed the switching. The 5-and 10-fold dilutions resulted in an increase in the molecular weight of the aggregates including the formation of fibrils both at 5 and 24 h, respectively (lanes 2-5; figure 5f ) as compared with the sample in 5 mM C12 FA (lane 1).

Ensemble kinetic simulation models validate the game-theoretic approach in elucidating
the dynamics of competing aggregation pathways

Parameter estimation
As mentioned in the Material and methods section, in the EKS model, four on-pathway rate constants (namely, k nu , k nu_ , k el , k el_ ), 10 off-pathway rate constants (namely, k con , k con_ , k nuf , k nuf_ , k el1f , k el1f_ , k el2f , k el2f_ , k fagf , k fagf_ ) and two off-on switching rate constants were considered (note that the forward and backward rate constants of switching each oligomer was considered the same, leading to only two switching parameters that need to be estimated, i.e. k swi , k swi_ ). Additionally, one needs to estimate two constants: p (which is simply a mapping constant that distinguishes the contributions of on-pathway oligomers from off-pathway oligomers to the ThT signal intensity) and pseudo-micelle concentration (concentration of the fatty acid near its CMC denoted by L). Following our published model in [34], the pseudo-micelle concentration was additionally estimated and not calculated directly from the FA concentration values at the CMC, since precise concentrations of pseudo-micelles are difficult to determine experimentally (only a fraction of total fatty acid concentration) as they are in dynamic equilibrium with other phases of micelle formation. This increased the number of parameters needed to be estimated to 18 from the EKS simulations. The potential complication is mitigated by the fact that our on-and off-pathway rate constants can be estimated separately using the respective control data. This makes it less cumbersome to estimate the remaining four rate constants (i.e. the two off-on switching rate constants k swi , k swi_ , the mapping constant p and the pseudo-micelle concentration L) from this offon switching dataset by significantly reducing the number of free parameters. A large parameter space from 10 −6 to 10 8 units with multiples of 10, was swept to estimate the value of each of the two switching royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191814 rate constants. Similarly, the pseudo-micelle concentration was varied from 0.01 to 1 units (with steps of 0.01), and p was varied from 10 5 to 10 8 units (with steps of 10 5 ). The estimated parameter values corresponding to the best fits are shown in table 2 in Appendix C. The benchmark on-and off-pathway rate constants (estimated separately from control data), were used to estimate the switching rate constants and obtain a global fit to the experimental ThT curves and monomer ratio values estimated from SEC measurements. The average R 2 of the off-to-on data is 0.974 and that of on-to-off data is 0.981.

Numerical results
The switching rate constants were sensitive specifically in the off-on dataset. The experimental data could not be fit in the absence of the switching rate constants and only a handful of switching rate constant combinations allowed an acceptable fit; the switching rate constants corresponding to the best fit to the experimental data are reported in table 2 in Appendix C. This directly proves the switching of off-pathway oligomers to on-pathway oligomers through the switching pathways due to the dilution of the system. The EKS simulations were conducted in the same way as the experimental set-up. For the off-to-on switching ( figure 6), first, combined off-and on-pathway  Figure 6. Correspondence between experimental results and EKS models on switching of pathways. (a) and (b) Experimental data (scatter data) on on-to-off and off-to-on pathways reproduced from figure 5b and f, respectively. Models based on EKS are shown as black lines. The intervention time points of 3 and 24 h (for (a)), and 5 and 24 h (for (b)) are shown as arrows. Panel (c) shows a phase diagram from EKS model at saturation, similar to figure 4 based on variations of the first two bridges. Here, the oligomer ratio of on-pathway to off-pathway was plotted as a heatmap (brighter colour, yellow, denotes on-pathway dominance while darker colour, blue, denotes off-pathway dominance) where the x-axis is bridge rate constant k con and y-axis is switching rate constant k swi . The phase diagram shows a dominance of on-pathway at low values of k swi and k con and dominance of offpathway for high values of k swi and k con .
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191814 simulations were executed, up to the switching time-point (of 5 or 24 h); all oligomer concentrations were noted until this point and they were then recalculated based on the amount of dilution at the switching time-point from the experiments. These altered concentration levels for each oligomer were next considered as the initial concentration of the combined off-and on-pathway simulation. Note that the second phase of the off-to-on dataset (figure 6a) did not show any lag time as can be seen in a usual unseeded on-pathway aggregation. Our model predicts a large conversion of offpathway species to on-pathway oligomers which results in a rapid formation of on-pathway fibrils (denoted by F ).
For the on-off dataset (figure 6b), stand-alone on-pathway simulations were executed exclusively up to the switching time-point (24 h) and the current oligomer concentrations were noted. These concentrations were then used to restart the combined on-and off-pathway simulation in addition to the pseudo-micelle concentration (that was also estimated in the parameter search step as an independent variable). Surprisingly, we found that the on-to-off-pathway dataset could be fit to our model both considering the switching rate constants, and also in the absence of switching rate constants generating comparable R 2 values; in other words, the switching rate constants had low sensitivity to the on-off experimental dataset. Probably as the on-pathway reactions are slow, very little on-pathway oligomers are formed at the switching time-point; as a consequence, this made the switching reaction flux slower than the previous case of off-on switching system resulting in overall lower sensitivity of the switching rate constants to the ThT data points from the experiments. While this precludes precise characterization of on-off switching, we do observe an overall decrease in fibril concentration compared with control data showing at least a qualitative impact of the switching reactions that convert the on-pathway oligomers into off-pathway species. Furthermore, we have also compared the phase diagram of the EKS model by plotting the oligomer ratio of on-pathway to off-pathway (as a heatmap) with varying bridge parameters and switching parameters during the saturation phase (75 h) (figure 6c). The total oligomer count scaled by their size from each pathway was used to compute this ratio. In this heatmap, brighter colour (yellow) denotes a dominance of on-pathway, whereas darker colour (blue) denotes a dominance of off-pathway. By doing so, four phases similar to those obtained from ROM were observed. For a low bridge and switching parameters, a dominance of on-pathway species was observed, whereas for a high value of bridge and switching parameters a prevalent off-pathway was observed (figure 6c); the light yellow and light blue regions depict the mixed pathway zones where both on-and off-pathway oligomers coexist. Note that a one-to-one correspondence between the phase diagrams generated from EKS and ROM models is not possible since the EKS models were built considering a detailed set of reactions, whereas the ROM models correspond to more bulk reactions involving fewer species.

Discussion
The data presented here is a first attempt in deciphering the complex phenomenon of protein aggregation pathways using a competition-based approach based on classical game theory. Aberrant protein aggregation is sensitive to environmental factors that determine the outcome of the aggregates [38,49]. Using the Aβ-fatty acid model system, we have employed a competition-based framework on simplified ROMs to gain preliminary insights. The results re-confirmed our previous observation that fatty acid concentrations modulate Aβ aggregation pathways [34]. Additionally, we discovered that the adoption of on-or off-pathway aggregates tightly depends on a set of rate constant ratios, which in turn suggest the thermodynamic stability (equilibrium constants) of the emerging aggregates. Moreover, α parameters are sensitive to the pseudo-micellar surfactant concentrations, L, which hold the key in modulating pathways. The models also provide insights into the feasibility of bridging pathways as a function of emerging higher-order aggregates. For example, the reduced order, sixspecies model predicts four different scenarios or dominant pathways of reactions which are strongly dependent upon the bridge, while also suggesting that α 4 is the key to the formation of larger aggregates in off-pathway. Stability arguments also show the larger aggregates in this system to be more stable (see Appendix B). The EKS simulations display a similar outcome; the simulations indicate that the larger the oligomer, the more significant the impact of that bridge upon formation of the respective fibril.
In our experiments, we note that the propensity to switch pathways is highest when the order of aggregate is the lowest (low-molecular weight) and increasingly becomes weaker as we move toward royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191814 higher-order aggregates along either pathway. This is in agreement with the theoretical studies noting the fact that in experiments, high molecular weight species refer to fibrils, while low-molecular weight aggregates then refer to the range of oligomers taken up in EKS and ROM simulations. Perhaps, a significant outcome of this study is the ability of the model to predict the emergence of oligomers by a set of kinetic and thermodynamic parameters, from a 'win' or 'lose' perspective (figure 4). Another key observation is the presence of multiple (neutrally) stable pathways in addition to simplistic onand off-pathways (figure 4; see Appendix B). The hybrid pathways, especially the off-on domain shown in yellow in figure 4, provide a range of possibilities for α 3 and α 4 to draw the aggregation dynamics away from toxicity. This is particularly significant for possible intervention strategies, pointing new lines of experimental and theoretical inquiries in the future. Numerical simulations and experiments both clearly support this qualitative result (figures 5 and 6), by revealing dilution to be a clear strategy to force the off-on transition.
Results of ROM indicate that when the ratio of non-bridge forward to backward rates is close to unity, the model achieves equilibrium quickly. However, when forward rates are considerably higher than backward rates, as is to be expected under normal circumstances, the system takes considerably longer to achieve equilibrium due to a cycle of over-shooting of species sizes resulting from a large difference in reaction rates. Similarly, corresponding EKS simulations indicate a 10 18 -fold difference in the forward and backward switching rate constants (table 2 in Appendix C) pointing to potentially irreversible effects of switching oligomers between pathways although the system may take a longer time to achieve equilibrium; this observation, however, pertains to our reaction system with fixed initial monomer concentration and is expected to show fluctuating dynamics by considering monomer or pseudo-micelle entry rates and stochastic effects of the switching of oligomers between the pathways.

Conclusion
The results presented here showcase the applicability of game theory on understanding amyloid aggregation pathways. This is significant because it provides an ability to predict the emergence of aggregates along multiple pathways along a temporal and equilibrium landscape map. Such a map can be further refined to see how it evolves as a function of a given interacting partner of Aβ, such as fatty acid as demonstrated here. A significant impact of this work could be realized with the potential for the prediction of the emergence of oligomers, which provides a handle for understanding the conditions at which toxic strains are formed and disappear. The simplified model presented here can be further fine-tuned into more sophisticated models by including more species along pathways, additional pathways and more interacting partners that can modulate the pathway, etc. In sum, the results presented here establish a new paradigm in understanding the complex dynamics of Aβ aggregation and provide impetus towards deciphering amyloid pathogenesis along with making therapeutic and diagnostic advances for such debilitating diseases in the future. ROM is, therefore, to be seen as a toy problem, which permits detailed analysis in a manner not possible with the very complex EKS model, where the parameters are all fixed (obtained through curve fitting with the experimental results). The results of the ROM provide insights into the system and help us ask the right kinds of questions about the kinds of experiments and EKS computations that needed to be performed. While some of the figures in the main text are restricted to the case n = 4, m = 20, these choices are made without loss of generality. Other choices of n and m (figure 7) have also been explored and the outcomes are seen to be qualitatively very similar to the one shown. Figure 7 shows the change in the % of the phase space (0 < ∝ 3 < 2, 0 < ∝ 4 < 2) taken up by each of the four pathways with changes to the pair (n, m). The bar graphs reveal these phases to barely change showing their ubiquity and theoretical significance. Appendix B, which focuses on the stability analysis of the system also reveals similar results.
Appendix B. Linear stability analysis of reduced-order model A linear stability analysis was conducted to confirm the conditions under which equilibria are stable and the sensitivity of these solutions to the parameters in this problem. We use the variables X 1 , Y 1 , X n , Y n , X m and Y m to represent the concentrations of the various perturbed species, while the equilibria for the monomers and oligomers from the two pathways are indicated by means of an 'e' in the subscript (i.e. B k,e represents the equilibrium concentration for the oligomer of size k). The central idea behind the stability analysis being that a stable equilibrium requires that the perturbed quantities eventually vanish, under certain conditions. The linearized system of equations for the perturbations yields the matrix, given by ðB 1Þ whose eigenvalues (denoted λ i , where i = 1-6) are indicative of the stability of the systems. Of key interest is the effect of the bridge parameters, α 3 and α 4 , and their effect on the stability of each model. A sampling of this effect is captured in figure 8, which depicts the contour plots of the eigenvalues of the Base Model for 0 ≤ α 3 , α 4 ≤ 2 for the special case when n = 8 and m = 40. In this figure, the lighter shades depict regions of low stability while the darker ones are more stable. The eigenvalue λ i corresponds to neutral stability. Overall, we find that the stability profile for equilibria corresponding to the Base Model does not change much for variations in values of n and m. The stability picture for the Base Model, however, is significantly different from that of Model 2. In the Base Model, one of λ 4 or λ 5 , is always zero for all values of α 3 and α 4 , while for Model 2, we observe switching behaviour between λ 2 and λ 4 , i.e. Model 2 shows greater sensitivity to the values of the bridge parameters. We distinguish two different characteristic effects, namely switching and crossing of the eigenvalues as the two bridge parameters are varied ( figure 9). The switching indicates a sudden, drastic change in behaviour of the species, where the course of domination of one species over other is abruptly reversed while the crossing is a more gradual version of this shift. In previous work [34], the switching has been compared with a sort of transcritical-like bifurcation in the system. Table 1 shows the switching and crossing points for eigenvalues λ 1 and λ 2 as α 3 and α 4 vary. As can be seen, there is crossing where α 3 = α 4 , whereas switching has an exponential relationship between the two parameters. A regression model indicates that switching occurs according to a 4 % 2:04Â1:55 a3 with R 2 = 0.953.
Simulations for other reaction rate regimes over a larger range of values for α 3 and α 4 (greater than 2) showed the switching and crossing to persist, as indicated by table 1. In previous studies with lowerdimensional models [34], we have seen such switching to occur as well, which appears to be indicative of the sensitivity of the system to the various pathways and species in the model and activation of one of these pathways under appropriate conditions. For the Base Model, with n = 2 and m = 4, we have λ 3 = 0 for low α 3 and α 4 , but as we increase these two parameters, λ 2 vanishes and then finally, with further royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 7: 191814 increase, λ 4 goes to zero. Thus, for any rate environment, the stability of the system near the point of equilibrium was found to be neutrally stable for sufficiently large values of the bridge parameters. The impact of species-size upon stability was also examined by studying the cases of (n, m) equal to (2,4) and (8,40) in addition to the standard case (4,20). Overall, in the case of the Base Model, no switching of eigenvalues was observed for a given species size environment as α 3 and α 4 were varied. For instance, in the Base Model, λ 5 is zero for all values of α 3 and α 4 when n = 2 and m = 4, and λ 4 is zero when n = 8 and m = 40. However, for Model 2, there is switching between λ 2 , λ 3 and λ 4 when n = 2 and m = 4, whereas when n = 8 and m = 40, we observe switching between λ 1 , λ 2 and λ 3 . In general, as n and m are increased, the overall magnitude of stability increases, i.e. the larger the species, the more stable the individual oligomer and also the overall system, appears to be.

Appendix C. The ensemble kinetic simulation model
Differential equations of the species: On-pathway species:  Figure 9. A sample case of eigenvalues depicting switching and crossing as a function of the bridge parameter α 3 . The eigenvalue λ 2 is denoted in blue while eigenvalue λ 1 is indicated in lavender.