Dynamical Systems Models for Plasma Dilution

Abstract

Introduction:

Recent experiments provide evidence that diluting the blood plasma restores the plasma environment to a more youthful level at least partially restoring the health of organs and tissues throughout the body.

Objectives:

We analyze whether a dynamical-systems model representing the plasma constituents could support the optimization process and help determine the appropriate dilution level, frequency and any simultaneous plasma infusions to achieve the most favorable outcome.

Methods:

We use a combination of a gradient descent, a simulated annealing and a genetic algorithm to find a population of models that fit illustrative data.

Results:

We analyze this population and present distributions of the model parameters and include a collection of plots of the dilution process for illustrative models.

Conclusions:

We find that, if we had more data to disambiguate the models, these dynmical-systems models could give useful predictions for plasma-dilution/addition experiments as well as support the interpretation of their results.

In a stunning pair of experiments, Mehdipour etal.1,2 showed that simply diluting the blood plasma with a combination of saline and 5% albumin, which they called neutral blood exchange (NBE), was sufficient to rejuvenate all three germ layers, improve cognition, and attenuate neuroinflammation in old mice. In fact, they found that a single NBE treatment, replacing 50% of the plasma in old mice, resulted in enhanced muscle repair, reduced liver adiposity and fibrosis, and increased hippocampal neurogenesis to at least the level seen in earlier heterochronic parabiosis experiments3 and surpassing the level seen in earlier heterochronic blood exchange experiments.4 Moreover, they showed that a single NBE treatment surpassed ABT 263, a senolytic that does not pass the blood–brain barrier, in its rejuvenative effect on the brain.

Additionally, the authors found that the protein composition of the plasma was reset to a more youthful, more rejuvenative composition, and although in this more youthful plasma, some proteins were reset to a lower level, as might be expected from a dilution experiment, perhaps surprisingly, many were elevated as a consequence of the NBE treatment. They noted that this was evidence that some of the age-elevated protein factors suppressed the age-reduced protein factors, and that once those suppressive proteins were diluted, the inhibited proteins were allowed to upregulate to a higher, more youthful level.

Noting that an equivalent procedure in humans, therapeutic plasma exchange (TPE), has already been approved for treatment by the US Food and Drug Administration (FDA), they further measured the blood proteomics of aged humans before and after TPE treatments and found that the protein composition reset in an analogous way to the old mice. The resetting and rejuvenating effect of TPE was further summarized in Refs. 5–9. Although TPE appears capable of inducing a resetting of the plasma composition, it is not yet fully clear how much plasma must be removed in order to see this resetting. The NBE experiments outlined above replaced approximately 50% of the plasma with saline and albumin, while typical TPE treatments replace the entire plasma volume.

Before continuing, we should also note that some researchers have considered the rejuvenating effect of some factors in young plasma and have performed some experiments where, rather than diluting plasma, they infuse young factors into the elderly.10–13 Of course, both may have merit, and the optimal protocol could very well involve a dilution combined with an infusion.

Beyond the measurement of the rejuvenative phenotype and the plasma proteomics, Mehdipour etal.1 also built a dynamical-system toy model that illustrated how a dilution might cause an enhancement to some plasma species in the system. This illustrative model had three species with the following set of differential equations (DEs): $$A\prime (t)=-\frac{1}{100}C(t),$$(1) $$B\prime (t)=-\frac{1}{20}A(t),$$(2) $$C\prime (t)=\frac{1}{10}A(t)-\frac{1}{100}C(t).$$(3)

We can see that in their toy model, species A positively affected C but negatively impacted B, while C negatively affected A and itself (which could be interpreted as its own decay). Schematically, this can be written as (4)

With this model, beginning with the initial values $$A(0)=200,B(0)=400,\text{and}C(0)=700,$$(5)

which they took to be diluted from $A_{BD}=1000,B_{BD}=400$ , and C_BD = 700, where the subscript BD stands for before dilution. Using this set of DEs and initial values, we reproduce their Fig.6 here in our Fig.1. As can be seen, this set of DEs gives interesting dynamics illustrating the relative values of some species decreasing (A and C), while another increases (B). It also displays oscillations on the way to equilibrium as the species overshoot before the feedback loops and decays remove the oscillatory energy, and it reaches equilibrium.

FIG. 1.

A reproduction of Fig.6 from Ref. 1 using the system of differential equations in Eqs. (1) through (3) and initial values given in Eq. (5).

Although the toy model achieved the illustrative purposes of Ref. 1, we would like to improve on their model in a few important way. We first note that their model only contained one steady state for A and C, and, therefore, the proposed initial state before the dilution with $A_{BD}=1000,B_{BD}=400$ , and C_BD = 700 does not correspond to a steady state of their model. In fact, a steady state for this system occurs when all the derivatives are zero, and the right-hand side of Eqs. (1) through (3) vanishes. This only happens when $A=C=0$ . On the other hand, B can take any steady state value and depends on the initial conditions. It “freezes out” when A vanishes. To improve on this, we would like to find models with at least two nontrivial steady states representing the before-dilution and after-dilution states of the system.

Another shortcoming of this toy model is that the species concentrations for A and C become negative before reaching the vanishing endpoint. This is clearly unphysical, and we would like to enforce that the model species never become negative in a more realistic system. In order to attempt an exploration of better models, we designed a system of DEs as described in Appendix A and wrote a computer code as described in Appendix B. We used our code to generate a population of models that fit the following data: $$\left(\begin{array}{c}{m}_0\\ m_1\\ m_2\\ m_3\end{array}\right)=\left(\begin{array}{c}0.273\\ 0.281\\ 0.552\\ 0.653\end{array}\right)\stackrel{\mathit{dil}}{\to }\left(\begin{array}{c}0.137\\ 0.141\\ 0.276\\ 0.327\end{array}\right)\stackrel{\mathit{evo}}{\to }\left(\begin{array}{c}0.477\\ 0.705\\ 0.300\\ 0.279\end{array}\right),$$(6)

where dil stands for dilution and evo stands for evolution. These values are meant to be illustrative with m ₀ and m ₁ increasing to a higher value after dilution and m ₂ and m ₃ decreasing to a lower value after dilution. We modeled these values on the proteins ROBO4, Angiogenin, VEGF-C, and Epiregulin from Fig.5E from Ref. 1. However, these values are only used in a demonstrative way in order to show what might be possible had we a robust dynamical model for plasma dilution. As we will see, we do not yet have enough data to disambiguate the many models that would fit the data.

The remainder of this article is organized in the following way. In the “Minimal Models” section, we describe the statistical properties of a population of minimal models which only contain the “measured” species in Eq. (6), and no extra species beyond these measured species. In the “Dilution Variations” section, we consider modified dilutions in the models described in the “Minimal Models” section. That is, we consider other dilution scenarios such as flat dilutions at different values, flat dilutions with a single or two species enhanced, and other possibilities. In the “Summary and Conclusions” section, we summarize our findings and conclude. We also include several appendices. In Appendix A, we briefly review dynamical models and describe the structure and parameters of the models used in our analysis. In Appendix B, we describe our algorithms for generating a population of models that fit the data. In Appendix C, we present plots illustrating the behavior of the minimal models following the standard dilution (the models and dilutions analyzed in the “Minimal Models” section). In Appendix D, we show plots illustrating the behavior of one of the minimal models under different modified dilution scenarios.

Minimal Models

In this section, we consider minimal models with only the measured species, but none of the extra species described in Appendix A. A selection of the models usedfor the analysis in this section can be found in Appendix C. Although the genetic algorithm attempted to add extra species (of the type u_i , h_i , and c_i —see Appendix B), none of the models with extra species had low enough model scores to remain in the population. However, it is, of course, possible for models with extra species to fit the data as well or better than these minimal models. In fact, in a separate calculation, we forced the model to retain at least four extra species and found many models that fit the data as well or better than the present minimal models. Nevertheless, we chose to present these minimal models because, first, they are minimal, and no extra species were required, and, second, since the models with extra species are more complex and have a much larger space of parameters to search, we were unable to obtain a large enough population of extra-species models with sufficiently low model scores that we could have confidence in the statistical properties of the distributions.

Each new minimal model took a little over 6 cpu-min to optimize, on average (some took much longer and some took much shorter). We generated a total of approximately 350 thousand new models that took a total of approximately 1,500 cpu-days of computation. We kept the 1,000 models with the best model scores (see Appendix B). In Fig.2, we show a distribution of the model scores in the final population. The maximum model score in the population depends on how many models we try. Greater run times produce models with lower model scores and lower population maximum score. In the future, with much greater data (especially data with more temporal measurements and data with a variety of dilution values), we would be most interested in only the few very best model scores and ignore all the rest. In the absence of sufficient data to disambiguate the many possible models shown here, we instead focus on traits shared by the many models that fit the data.

FIG. 2.

A distribution of model scores from our best 1,000 models with only four measured species.

As described in Appendix B, for each species’ equations, when creating a new model from scratch, we randomly chose (with a flat distribution) between 2 and 7 terms [the number of terms on the right of Eq. (A6)], including the final decay term. However, the genetic algorithm sometimes took a preexisting model from the population and either added a term or removed a term (see Appendix B), allowing the number of terms to exceed these bounds and change the shape of the distribution. Therefore, the final shape of the distributions is a combination of the flat distribution convoluted with the mating, the mutations, and the survival of the fittest algorithms. With this in mind, we present a distribution for the number of terms in the equation for each constituent in the top plot of Fig.3. We can see that the number of terms per equation is peaked at two terms for m ₀ and m ₂, at three to four terms for m ₁, and at four to five terms for m ₃. As expected, there were no models with only one term (only the decay term) on the right of any equation, since the steady state for such a species would necessarily be zero. There must be something to balance the decay in order for the species to be nontrivial. After the dilution, m ₂ rises a very small amount, which apparently only requires one additional term, m ₀ increases a moderate amount, which also only requires one additional term, but m ₁ increases twice as much as m ₀ and apparently needs at least two or three additional terms. On the other hand, m ₃ decreases following dilution, the opposite of the other three, and this requires three to four additional terms. Furthermore, as we can see in the middle of Fig.3, typically, only the decay term is negative for m ₀ and m ₂. Once again, this must be the case, else m ₀ and m ₂ would be trivial with both terms negative and no nontrivial steady state. On the other hand, one or more of the extra terms is typically negative for m ₁ and m ₃, presumably due to the extra complexities of their evolutions. For each species’ equation, we have also analyzed how many of their terms are of the form Eq. (A2), with only a single species, or of the form Eq. (A4), with two species, and show a plot in the bottom of Fig.3. Once again, m ₀ and m ₂ are similar with around 80% being single-species terms, while m ₁ has approximately half of its terms with two species and approximately 80% of the terms for m ₃ have two species.

FIG. 3.

These plots show the fraction of models or terms that have certain properties in the equations for each species. The blue, orange, green, and red lines represent the equations for $m_0,m_1,m_2$ , and m₃, respectively. The top plot is for the number of terms per equation. The middle plot is for the number of negative terms per equation. The bottom plot is for the number of species per nondecay term for each equation. In other words, whether a term was of the form Eq. (A2) or Eq. (A4).

If we look at each equation in greater detail, we begin to see further patterns. In Fig.4, we plot which species is included for single-species terms for each equation. These are terms of the form Eq. (A2). The blue line is for the first equation that gives ${\dot{m}}_0$ . We see that a little over 80% of its single-species terms contain m ₁, and a little more than 10% contain itself. The orange line, on the other hand, represents the equation for ${\dot{m}}_1$ , and nearly 70% of its single-species terms contain m ₀, while nearly 30% contain itself. It makes sense that these two are tightly coupled since they share the property that both increase following dilution to a steady state that is higher than the previous steady state. As we can see, species m ₂ and m ₃ are similarly related for single-species terms. We see from the green line that approximately 90% of the single-species terms for ${\dot{m}}_2$ contain m ₃, while approximately 80% of the single-species terms for ${\dot{m}}_3$ contain m ₂, in red.

FIG. 4.

Fraction of single-species terms containing a particular species. The blue, orange, green, and red lines represent ${\dot{m}}_0,{\dot{m}}_1,{\dot{m}}_2$ , and ${\dot{m}}_3$ , respectively. The horizontal axis represents whether the species in the single-species term is $m_0,m_1,m_2$ , or m ₃.

In Fig.5, we plot distributions of the coefficient and thresholds for the most important single-species terms in each equation. These lines correspond with the peaks in Fig.4. Once again, ${\dot{m}}_0$ and ${\dot{m}}_1$ have similar properties. For ${\dot{m}}_0$ , the coefficient takes the value $g_{01}=1.00\pm 0.42$ , and the threshold takes the value $K_{01}=0.76\pm 0.25$ . Meanwhile, ${\dot{m}}_1$ has $g_{10}=1.05\pm 0.54$ and $K_{10}=0.714\pm 0.44$ . The main difference is that the peaks for g ₀₁ and K ₀₁ are narrower and taller, especially in the case of the threshold. The thresholds are notable for their proximity to the final stead-state values for m ₀ and especially m ₁, so that the steady state occurs when they are near their steady state. Moreover, the threshold effect is important for these species as they traverse a large fraction of the transition from low species concentrations to half of their saturation value. The coefficient for ${\dot{m}}_2$ is similar to ${\dot{m}}_0$ and ${\dot{m}}_1$ with $g_{23}=1.07\pm 0.48$ ; however, ${\dot{m}}_3$ stands out sharply from the others with $g_{32}=0.69\pm 0.44$ . Presumably, this is because m ₃ needs to decrease following dilution and, therefore, needs to reduce this positive feedback since it would tend to cause m ₃ to increase. However, it should be remembered from Fig.3 that single-species terms are only a small contribution to ${\dot{m}}_3$ . Furthermore, their thresholds are $K_{23}=1.12\pm 0.36$ and $K_{32}=1.07\pm 0.57$ , which are very similar, although K ₂₃ is narrower. Interestingly, the data appear to favor a threshold for m ₂ and m ₃ that are well above both their initial steady state and their final steady state. This suggests that the threshold may not be very important, and they spend the majority of their time well below saturation. Therefore, thepositive feedback from these terms is nearly linear.This is interesting since both have relatively little change between their diluted value and their final steady state.

FIG. 5.

Percent of single-species terms with coefficients (top) and thresholds (bottom) with given values (horizontal axes). The coefficients and thresholds refer to g and K, respectively, in Eq. (A2). For each equation ( ${\dot{m}}_0,{\dot{m}}_1,{\dot{m}}_2$ , or ${\dot{m}}_3$ ), we have only shown the coefficients and thresholds for the term that occurs most frequently. That is, for the peaks in Fig.4.

In Fig.7, we plot which species are included for double-species terms for each equation. These are terms of the form Eq. (A4). The first thing that we note is that, in the equation for ${\dot{m}}_j$ , where j = 0, 1, 2, or 3, the majority of double-species terms contain themselves. In other words, they contain a product of the form $m_j{m}_k$ , where $k\ne j$ . So, for example, approximately $71\%,96\%,70\%$ , and 68% of the double-species terms for ${\dot{m}}_0,{\dot{m}}_1,{\dot{m}}_2$ , and ${\dot{m}}_3$ contain a product of the form $m_0{m}_k,m_1{m}_k,m_2{m}_k$ , and $m_3{m}_k$ , respectively, where k is some other species. For the equation for ${\dot{m}}_0$ , the other species in $m_0{m}_k$ is spread nearly evenly with between 20% and 30% for each of m ₁, m ₂, and m ₃. For ${\dot{m}}_1$ , on the other hand, a dominant 60% of the double-species terms contain $m_1{m}_3$ , where m ₃ is the only species that decreases following dilution, while m ₁ is the species that increases by the most. Additionally, 34% of the terms for ${\dot{m}}_1$ contain $m_1{m}_2$ , where m ₂ is the other species whose final steady state is lower than the initial one. It appears to be important for m ₁, which grows substantially, to have a negative (see the next paragraph) feedback mechanism linked to the species that decreases. On the other hand, nearly no terms have $m_1{m}_0$ , where m ₀ increases following dilution. Turning to ${\dot{m}}_2$ , 40% of the double-species terms contain $m_2{m}_3$ , perhaps because they share the property that their final steady state is lower than their initial one. Finally, for ${\dot{m}}_3$ , its most frequent double-species term is $m_0{m}_1$ , which interestingly does not contain m ₃. Apparently, the decreasing m ₃ needs a positive (see the next paragraph) feedback depending on the levels of the increasing m ₀ and m ₁. However, there are nearly as many equations with a term of the form $m_2{m}_3$ and $m_1{m}_3$ , both of which do contain m ₃. This unique feature must be due to the unique evolution of m ₃ following dilution.

In Fig.6, we plot the coefficients of the most important double-species terms. As can be seen, we do not include any double-species plots for ${\dot{m}}_0$ or ${\dot{m}}_2$ . The reason is that the statistics are quite poor as can be deduced from Fig.3 combined with Fig.7, and, therefore, these terms are not present in a large fraction of models. It appears that ${\dot{m}}_0$ and ${\dot{m}}_2$ rely mainly on single-species terms. These are the two species with the smallest increase following dilution. From these same figures, we can also see that ${\dot{m}}_1$ uses double-species terms much more frequently, and that half of its terms contains a product of two species. Presumably, this is because it increases the most drastically following dilution and needs to be closely coupled with other species in order to time both its increase and its approach to equilibrium. ${\dot{m}}_3$ uses, by far, the greatest number of double-species terms at greater than 80% of its terms. This must be due to the fact that it is the only term that decreases following dilution. Again, we think that these terms better enable it to time its evolution to a new steady state. For these reasons, we focus on ${\dot{m}}_1$ and ${\dot{m}}_3$ .

FIG. 6.

Fraction of double-species terms with coefficients of given values (horizontal axis). The thresholds for these terms are shown in Fig.8. As described in the text, we have only shown the coefficients for the terms that occur most frequently.

FIG. 7.

Fraction of double-species terms containing two species. The plots are, from top to bottom, for ${\dot{m}}_0,{\dot{m}}_1,{\dot{m}}_2$ , and ${\dot{m}}_3$ . The horizontal and vertical axes are for the two species in the product.

As we look at the (blue and orange) distributions for ${\dot{m}}_1$ in Fig.6, we see that both are negative. In particular, $g_{112}=-0.91\pm 0.79$ and $g_{113}=-1.05\pm 0.65$ . This makes sense since, as we can see in Fig.5, the single-species terms tend to feedback positively, causing it to increase its amount. In order for this positive increase to end, ${\dot{m}}_1$ needs a negative feedback mechanism that keeps it from increasing too much. This negative feedback should be sensitive to both its own value and the value of other important species, especially to m ₃. The coefficients for ${\dot{m}}_3$ (the green, red, and purple lines in Fig.6), on the other hand, are split between positive andnegative feedbacks. Namely, $g_{301}=0.76\pm 0.39,g_{313}=-1.14\pm 0.60$ , and $g_{323}=1.09\pm 0.65$ . For ${\dot{m}}_3$ , the desired effect is just the opposite. As we can see in Fig.5, the single-species terms tend very strongly to be positive feedback and, thus, tend to increase the concentration of the species. For this reason, it appears that ${\dot{m}}_3$ does not have an abundance of single-species terms. In fact, more than 80% of its terms are double-species terms, as noted previously. We will come back to this point in the next paragraph. In order to decrease following dilution, ${\dot{m}}_3$ relies on joint negative feedback fromitself (as is necessary) as well as m ₁ (the red line). Since this term tends to decrease the amount of m ₃ along with its decay, it now needs positive feedback to compensate and allow it to reach a steady state. It does this with the combination $m_0{m}_1$ , which does not depend on itself and, therefore, must be positive as well as from $m_2{m}_3$ .

In Ref. 1, the authors noted that it was interesting that many constituents of the plasma increased rather decreased following dilution. This might be surprising since, naively, one might expect a dilution to reduce the concentrations to a lower steady state. Moreover, inspecting their Fig.5 reveals that all of the reported mouse proteins increased beyond their original values, while nearly 80% of the reported human proteins increased. This appears to suggest that it is more likely for a protein to increase following a dilution event rather than decrease, although they only showed a fraction of the plasma proteins, and it is possible that there is a bias in which proteins were shown. In this paper, we have shown that, at least within the constraints of the limited data we have, single-species terms tend to increase species concentrations, while a decrease tends to require greater contributions from double-species terms. We might hypothesize that the reason a decrease is less likely is that the evolutionary path to enzymes that were catalyzed by a single species was more efficiently explored, and that enzymes requiring two species for catalysis were more difficult to achieve stochastically and, therefore, less efficiently explored. Further analysis of dynamical models might shed greater light on this biological trait.

We note one last property of these coefficients before we move on. We can see that there is a preference for larger magnitude coefficients, around 1, rather than smaller ones. Although the percentage of terms with coefficients larger than 1 decreases rapidly, it should be kept in mind that every search was begun with a random value whose magnitude was between 0 and 1. The simulated anneal process and gradient method were free to search outside this range, but their initial starting point effectively reduced the distance they searched beyond unity. We attempted to increase the initial value to a magnitude that was randomly chosen between 0 and 2. What we found was a similar deficit around 0, followed by a quick growth to a maximum value around 1, a nearly flat distribution between magnitudes of 1 and 2, followed by a quick decline afterward. There are some exceptions, such as the green line for ${\dot{m}}_3$ in Fig.6 or the red line, again for ${\dot{m}}_3$ , in Fig.5, where the peaks occur at a value lower than 1, but generally speaking, this is the behavior. These same comments apply to the coefficients of the single-species terms and the decays that were at their peak near 1 which we turn to later in this section. We suggest this to be a property of the fact that each of these models represents an infinite class of models that differ by a scaling of their coefficients and decay rates. In order to keep our analysis simple and to focus on the classes of models rather than being overwhelmed by scaled examples within a class, we chose to keep the original random range for the magnitude of our coefficients and decay rates between 0 and 1.

In Fig.8, we show distributions of threshold values for the most frequent double-species terms. These are the Ks in terms of the form Eq. (A4). The main thing to note here is whether the thresholds are very small, intermediate, or very large compared to the concentrations of the species during the dynamic response to dilution. We can see that there are two instances of very small threshold values, and both occur for ${\dot{m}}_1$ and are both thresholds for m ₁ (top two plots). Specifically, $t_{1121}=0.57\pm 0.49$ and $t_{1131}=0.48+-0.40$ (the index stands for ${\dot{m}}_1$ with product $m_1{m}_2$ or $m_1{m}_3$ and the threshold for species m ₁) both with peaks far to the left of this and long tails to the right (the mean is deceptive in this case). When this threshold is very small compared to the concentrations, it signifies that the catalyst is nearly always saturated in this species, and that the reaction rate is nearly constant in this species (although not necessarily in the other species of the product). We remember from Fig.6 that the coefficients for these terms are negative, so that this represents the rate of destruction of m ₁. This probably could have been modeled as a single-species term without m ₁in the destruction term; however, our policy (see Appendix A) does not allow destructive terms unless the species being destroyed is part of the term. The reason for this was that a species has to be present in the catalyst in order for it to be destroyed, and, therefore, with this constraint, the species concentration cannot become negative. The threshold for the other species in these terms has a roughly flat distribution between intermediate and largevalues, namely, $t_{3122}=0.91\pm 0.51$ and $t_{3133}=0.90\pm 0.48$ , suggesting that saturation is not common. The dependence on the other species in this term is sensitively dependent on the concentration of the other species in some models, while it is nearly linear in the other species in other models.

FIG. 8.

Fraction of double-species terms with thresholds of given values (horizontal axis). The coefficients for these terms are shown in Fig.6. As described in the text, we have only shown the thresholds for the terms that occur most frequently.

Turning to ${\dot{m}}_3$ , we see that only one term has a strong preference for intermediate thresholds (the fourth plot of Fig.8). It is interesting that this is the only term with anegative coefficient (see Fig.6) and depends on bothitself m ₃ and m ₁. These thresholds are $t_{3131}=0.68\pm 0.41$ and $t_{3133}=0.65\pm 0.36$ , again with a peak to the left and a long tail to the right. As we have seen in previous paragraphs, these are typically the species that stand out, and this is no exception. This is the only major catalyzed destructive reaction for m ₃, and since the threshold value is intermediate for both m ₃ and m ₁, the rate of destruction varies greatly by the concentrations. In particular, directly following dilution, it is at its greatest rate and nearly linear. However, as m ₁ increases toward its final steady-state value, it approaches and then surpasses its threshold value, slowing down the rate until it is nearly constant in m ₁. This likely helps m ₃ to stabilize at its final value and not decrease too much. The dependence of the rate on itself m ₃, on the other hand, is nearly constant before the dilution when m ₃ is high but approaches linear in m ₃ after the dilution and remains so to the final steady state. This makes sense since, if it increased from its steady state value, it would want to pull back at a level commensurate to its distance from the steady state. The other thresholds appear to disfavor very low values but only have slight to no favor between intermediate and larger values, namely, $t_{3010}=0.88\pm 0.51,t_{3011}=0.96\pm 0.54,t_{3232}=0.86\pm 0.49$ , and $t_{3233}=1.01\pm 0.53$ . This suggests that a rate that is either significantly varying or nearly constant is adequate in the concentrations but a linear rate would be too much.

In Fig.9, we consider the decay rates for each species. The λ refers to the coefficient in Eq. (A1). The red line for m ₃ immediately stands out. It has a very strong preference for very small decay rates, ${\lambda }_3=0.14\pm 0.12$ . This suggests that its decrease following dilution depends on the species having a relative stability compared to the other species. The destruction of m ₃ apparently depends heavily on catalysts, and it does not tend to break apart on its own (on the scale of these experiments). The next lowest decay rate is for m ₁, ${\lambda }_1=0.46\pm 0.27$ , another term which has a strong destructive contribution from catalysts, and this is correlated with relative stability on its own. Although in this case, its uncatalyzed decay does happen at an appreciable rate on the time scale of the reactions, m ₂ and m ₀ have decay rates that are a bit higher, ${\lambda }_2=0.75\pm 0.34$ and ${\lambda }_0=1.02\pm 0.43$ , respectively. This is noteworthy because m ₂ and m ₀ are the two species that do not have strong catalyzed destructions. So, they make up for this by falling apart more easily, thus stabilizing the catalyzed creations of them in order to achieve a steady state.

FIG. 9.

Fraction distributions of the decay rates for each species.

In this section, we consider statistical distributions of models that fit the illustrative data in Eq. (6). Since the data for these four species were so limited, there was a very large amount of ambiguity in which model could fit the data. However, we have shown that some correlations could, nevertheless, be inferred that might be useful as more data are collected, and the models are improved. We also note that these models and the associated statistical inferences are further limited by the ability of the algorithms to generate instances of unique classes of models that fit the data. We have attempted to make a robust algorithmic system that does not prefer any type of model, but we can never be sure, and further studies are always warranted. Finally, due to our limited computational resources, we were necessarily constrained to including a very small set of measured species, only 4 in this analysis, and no extra species that were not measured. This is important because the presence of feedback interactions with other species not included here could modify the results of this study. It may be that the dependence of one species on another found here might be shifted to a dependence on other species which were not modeled here. This is a major direction we would like to consider in future studies by enlarging our computational power and analyzing models with a greater number of species. Moreover, we would like to include further data, both with a variety of dilution values and data at more time points. This last point, of course, will depend on the experiments being performed but will lead to the greatest improvement in these models.

Dilution Variations

In this section, we consider variations on the dilution from the experiment. We would like to explore what predictions a dynamical system trained on sufficient data might make. Since we do not have much data for training, we cannot expect the predictions of individual models to match Nature yet. Instead, we have a population of models and, therefore, a population of predictions, each unique. On the other hand, we can see patterns among the models, suggesting that there might be a greater likelihood of a prediction to have certain characteristics. Nevertheless, any predictions made at this point are premature and should be seen as illustrating what might be achieved had we sufficient data to find a more robust, realistic model of an actual biological system. We note that all the figures in this section have dilutions that range from 0.02 up to 1.4. Dilutions above 1.0 are, of course, not dilutions. They are enhancements (infusions). However, for convenience of language, when referring to them together with the dilutions, we will simply call them all dilutions. On the other hand, when we focus on values larger than 1, we will use the more appropriate term enhancements.

We begin by considering a flat dilution where we dilute every species by the same amount, although not necessarily the 50% discussed in the “Minimal Models” section. In Fig.10, we show the range of final values for each of the species if they are all diluted by the same amount, shown along the horizontal axis. The solid curves represent the median values, while the shaded regions show the range of the middle 67% of the models. We can see that if we dilute all of the species by any value below 1, 67% of the models evolve to the same final steady state that occurs with a 50% dilution. We will call this the original final steady state (OFS). This suggests that in these simple models, the OFS is an attractor, and we will see further evidence of this in the rest of this section. If we instead enhance the species, we see that the majority of the models tend toward final values where m ₀ and m ₁ are nearly 0, while m ₂ tends toward a range between 0.98 and 1.4 and m ₃ tends toward a range between 2.3 and 11.3, all in 67% of our models.

FIG. 10.

Final values following a flat dilution. The shaded region covers the final value of 67% of the models. All the species were diluted by the same amount, shown along the horizontal axis. The dotted lines give the final steady states of the original data (the OFS).

In Fig.11, we consider a flat 50% dilution for all the species except one. The excluded species is instead diluted at a range of values including greater dilution, lesser dilution, and enhancement. If the dilution or enhancement is greater than 50% for the decoupled species, this could be achieved by a dilution followed by an infusion of the separated species. However, if the singled-out species is diluted by a greater amount than the rest of the plasma, this could not be done with a dilution–infusion combination. Rather, this might require the infusion of some agent that could bind to the species and remove it in some way. These comments also apply for other scenarios discussed in this section.

FIG. 11.

Final values following a 50% dilution of all species except for one. The single species not diluted by 50% had a varying dilution given by the horizontal axis. The shaded region covers the final value of 67% of the models. The dotted lines give the final steady states of the original data (the OFS).

We begin with the top plot, which considers a decoupled m ₀ dilution. We see that, given a 50% dilution for everything else, the system is very insensitive to the dilution value for m ₀. For nearly the entire range, the final value coincides with the OFS for 67% of the models. If the dilution reaches down below approximately 0.2, the final steady states begin to spread out, although the median value still coincides with the OFS. A very similar statement can be made for dilutions of m ₁ and m ₂. For the majority of dilution values, the final steady coincides with the OFS. However, if the dilution of m ₁ falls below approximately 0.3 or the dilution value of m ₂ falls above approximately 0.9, the final steady state begins to spread out, once again with the median remaining true to the OFS. These are the three species which increase following dilution. The main difference between m ₀ and m ₁ on the one hand and m ₂ on the other hand is that they are sensitive to opposite ends of the dilution spectrum. The reason for this is that the OFS for m ₀ and m ₁ is above their initial steady state; therefore, they are more sensitive to severe dilutions which are in the opposite direction of their OFS. The behavior of m ₂ is the opposite with its OFS below its initial steady state; therefore, it is more sensitive to enhancements, which lead in the opposite direction of its OFS.

The species that really stands out here is m ₃, which is sensitive to both dilutions above 0.7 and enhancements. Or to put it differently, in order to reach the final steady state of the data, m ₃ must be diluted by at least 30%. If it is above that, even with all other species diluted by 50%, the model will tend toward a different final steady state. m ₃ is the most unusual as we have seen throughout this note as it not only has a lower final steady state but also decreases following dilution and requires the greatest complexity of its dynamical equations. It appears that this confers a greater sensitivity on it than the other species. Interestingly, with enhancements, 67% of the models tend toward final values with m ₀ and m ₁ at or very near 0, while m ₂ and m ₃ tend toward values in the range 0.97 and 1.4 and 2.3 and 11.6, respectively. These are the same final values as occurred in Fig.10, where all species were enhanced by the same amount. This suggests that m ₃ might have been the dominant change causing evolution toward this final state. Having the other species diluted at 50% here helped the transition occur, so that m ₃ could take a larger range of dilutions and enhancements of 0.7 and above to achieve these final states compared with the flat dilutions shown in Fig.10.

We also considered manipulation of two species independently while keeping the remaining two at a 50% dilution. So, for example, we diluted m ₂ and m ₃ at 50% and then varied the dilution or enhancement of m ₀ and m ₁, independently, over the range 0.05–1.4. We found roughly what we would expect if we considered the combined effect of two modified dilutions that are shown in Fig.11. Since there are too many plots to display here and since there is a large degree of similarity between the plots, we will focus on a few illustrative examples here. In Fig.12, we present the median final value of m ₁. Along the horizontal axes and vertical axes, we present the dilution values for two species, while we hold the dilution of the other two species at a constant 50% dilution. So, although the final value after dilution is always presented for m ₁, the dilution values depend on the values of the axes or are 50% if not on the axes. In these plots, the yellow region shows where the median final value of m ₁ is the OFS, the right side of Eq. (6). The blue-green area is the transition region between final states, and the dark blue region shows the region where the alternate final steady state (AFS) occurs following dilution, shown on the right side of the bottom plot of Fig.11, for example. Although we only show plots for the final state of m ₁, the plots for the other species, given the same dilutions, are analogous. The exact positions of the lines separating the regions are slightly altered but is in the same region and has the same qualitative shape.

FIG. 12.

Distributions of median final value of m₁ when two species were diluted or enhanced, independently. These two species are specified in the horizontal and vertical axes. The other two species, not on the axes, are diluted at 50% in each case. FS stands for final steady state.

In the case that we dilute m ₂ and m ₃ at a constant 50% but vary the dilution of m ₀ and m ₁ (top plot of Fig.12), we find that, for the majority of dilution and enhancement values, the final state is the OFS. However, if we dilute either m ₀ or m ₁ to very small dilutions, the final state switches to the AFS. This is shown as the darker region at the bottom left. This is because, as we can see in the second plot of Fig.11, the final state begins switching at very low dilutions of m ₀ and m ₁, and when we combine very low dilutions of both, the switch becomes complete. We note that the solid orange line in the second plot of Fig.11 corresponds with a horizontal line in this plot at a height of 0.5 and is yellow for the entire length, as expected. Similarly, the solid orange line in the first plot of Fig.11 corresponds with a vertical line in this plot at a value of 0.5 and is also yellow for the entire height, as expected. Similar comments apply for the horizontal and vertical dashed lines in the other plots of Fig.12. We also note that we are only plotting the median value and do not have a clear way to include the 67% confidence region. In the second plot of Fig.12, we show the effect of varying the dilution of m ₀ and m ₂. We can see that the AFS switches from the bottom left to the bottom right. This makes sense since, as we can see in the third plot of Fig.11, the final state switches for high values of m ₂, which is the right side in this plot. Again, we find that low values of m ₀ hasten the transition. If we vary the dilutions of m ₁ and m ₂ instead (not shown), we find a qualitatively very similar plot to the second plot of Fig.12. Again, the transition region is at the bottom right. However, the transition line has a slope that is shallower.

The shape changes significantly when we vary the dilution of m ₃, as we would expect based on the bottom plot of Fig.11. If we also vary m ₀, we get the behavior shown in the third plot of Fig.12. The transition is most sensitive to variations in the dilution of m ₃. However, a smaller dilution value for m ₀ makes it more likely, while a larger dilution value (or enhancement) of m ₀ makes the transition less likely, causing the transition line to slope to the right. Varying the dilution of m ₁ along with m ₃ (not shown) creates a qualitatively similar plot with the slope of the transition line to the right, but the slope is shallower. Finally, if we vary the dilution of m ₂ along with m ₃, as expected, the slope of the transition line flips to the left. We expect this, based on the third plot of Fig.11, since the transition is sensitive to large values of m ₂ and, therefore, should occur more quickly the higher the dilution value of m ₂.

In Fig.13, we consider what would happen if we only diluted (or enhanced) one species but left all the others at their original value. We can see that the behavior splits into two groups. In the first, as shown in the first two plots, if we enhance either m ₀ or m ₁ alone, the model evolves toward the OFS. This is likely because the OFS is higher than the initial steady state for m ₀ and m ₁, and so this begins the evolution in the direction of the OFS. On the other hand, m ₂ and m ₃ are just the reverse. Their OFS is below the initial steady state; therefore, a dilution moves them in the direction of the OFS. As we see in the bottom two plots, any dilution of either of these causes an evolution to the OFS. On the other hand, if we perturb any single constituent of the group in the direction opposite the OFS, it evolves toward a different final state. Interestingly, that final state is essentially the same for any single-species perturbation and is, in fact, the same as when we enhance all four species (as shown in Fig.10) or when we dilute all species at 50% except for m ₃, which is enhanced (as seen in the bottom of Fig.11). Namely, in all four cases here, both m ₀ and m ₁ tend toward very small or zero values, while m ₂ tends toward a range of 1.0 and 1.4 and m ₃ tends toward a range of 1.8 and 10.5 in the first plot, 2.2 and 11.5 in the second, 2.2 and 10.4 in the third, and 2.3 and 11.4 in the fourth. Although the range for m ₃ is not exactly the same, the regions overlap significantly. This AFS region is interesting because it was not included in the training data, and we did not build it into the models. It appears that, within the limitations of the data that we used, the presence of the steady states given in Eq. (6) and the evolution between the two is sufficient to predispose the model to have another steady state in this region. At least, it occurred in 67% of the models obtained using our algorithm. That is, the other 33% of our models that fit the data do not have a steady state in this region. So, if we had more data, we might continue to find a steady state in this region or we might not. However, a more robust model, generated with a greater quantity of quality data, might predict the presence of other steady states that might be interesting and might suggest ways of achieving an evolution to the alternate steady state (ASS).

FIG. 13.

Final values following a dilution of only one species while not diluting the others. The single species diluted is given by the horizontal axis. The shaded region covers the final value of 67% of the models. The dotted lines give the final steady states of the original data (the OFS).

Another important feature of Fig.13 is that the original steady state is not stable. A perturbation in any direction of any single species is sufficient to drive the model away from the original steady state toward either the OFS or toward some other final state, as described in the previous paragraph. This property is likely not physical as we expect the levels of constituents in our plasma to be stable to small changes. This is a feature of our models we would like to improve in the future. However, we would prefer to not build in a region of stability by hand since we do not know the sizes of these regions physically, and they likely are different for each species being perturbed. Furthermore, we suspect that the number of steady states as well as their stability is impacted by the size of the dynamical model. That is, as we increase the number of species in our model, including unmeasured, hidden and conserved species, as described in Appendix A, the number of steady states will likely increase, and their stability structure could change. Moreover, we would prefer to let the data determine the features of the model. Given sufficient data of the right type, the steady states should take on more realistic properties. In particular, as a place to start, we would like to see dilution experiments with a range of dilutions, including dilutions where the plasma returns to its original steady state. The OFS, as given on the right of Eq. (6), on the other hand, is stable to a wide range of perturbations of a single species in our models.

In Fig.14, we show the result of diluting two species independently while holding the others at their original value. As in Fig.12, we have only plotted the median final value of the models with each dilution. Moreover, as in that previous figure, we see that the dilution space is divided into three regions: the OFS (in yellow on the left and dark blue on the right) and the AFS (in dark blue on the left and green on the right) and an intermediate transition region separating them. However, this time, there is a qualitative difference between the behavior of m ₀ and m ₁ on the one hand and m ₂ and m ₃ on the other hand. In particular, the transition region is broadened for the latter two. For this reason, we have plotted the distribution for m ₁ in the left column, and these plots are also representative of m ₀, while we have plotted the distribution of m ₂ (also representative of m ₃) in the right column. We have also included dashed lines when the dilution of either constituent is 1.0, representing no dilution at all. These dashed lines correspond with the solid lines in Fig.13. As expected, these dashed lines cross in the transition region, where none of the species are diluted, and the final state is the original starting state [the left side of Eq. (6)], which is in between the OFS and the AFS.

FIG. 14.

Distributions of median final value of m ₁ on the left and m ₂ on the right when two species were diluted or enhanced, independently. These two species are specified in the horizontal and vertical axes. The other two species, not on the axes, are held at their original value in each case. FS stands for final steady state.

In the top row of Fig.14, we see the behavior of the constituents when only m ₀ and m ₁ are diluted. As we would expect, based on first and second plots of Fig.13, the OFS is obtained if both species are increased, but the AFS is preferred if they are both diluted. Additionally, in the top left and bottom right regions, we see that an enhancement of one but a dilution of the other results in the OFS if the dilution is not too much, but transitions to the AFS if the dilution is large. We also see, in the right plot for m ₂ and m ₃, that the transition region takes up a larger fraction of the space in the top left and bottom right, when one species is enhanced and the other diluted. This is presumably because m ₂ and m ₃ are more sensitive to the level of dilutions/enhancements, which appears to be a consequence of their decrease following dilution. In the second row, we show plots for dilutions in m ₀ and m ₂. Once again, if we compare this behavior with the first and third plots of Fig.13, we see that enhancing m ₀ and diluting m ₂ result in the OFS, while the opposite results in the AFS, as expected. On the other hand, either enhancing both or diluting both results in one or the other depending on which is stronger, with a transition region separating them. Once again, the transition region is larger for m ₂ and m ₃ (in the right plot). The plots (and behavior) are very similar for dilutions of m ₁ and m ₂, and so we do not include plots for this combination here. However, we do note that the transition region is closer to the diagonal for the $m_1-m_2$ combination. In the third row, we consider the combined dilutions of m ₀ and m ₃. In this case, the dilutions and enhancements of m ₃ are overpowering, and the transition region is nearly a vertical line at a dilution of 1.0 for m ₃. We have seen throughout this note that the model is very sensitive to m ₃, and that this appears to be connected to the fact that it decreases not only following dilution but also from its diluted value, at least with the original dilution values. Nevertheless, we do see a little slant to the transition line, showing that enhancements of m ₀ do increase the OFS region a little, and dilutions of m ₀ decrease it a bit, as expected. The combined $m_1-m_3$ dilution behavior is very similar to this, and we do not show it. On the last row, we consider dilutions of m ₂ and m ₃. Since dilutions and enhancements tend to have similar effects, we see that dilutions of both give the OFS, while enhancements of both give the AFS. Enhancement of one and dilution of the other is split, depending on which is stronger. This case is the opposite of the first row, where we consider the $m_0-m_1$ case.

Summary and Conclusions

In this note, we considered dynamical systems models for plasma dilutions. The plasma dilution experiments, so far, have only been done at one dilution value, and the concentrations of plasma proteins have only been measured once before and once after the dilution event. Because these data are so limited, we were not achievable a unique model that faithfully represents the behavior of the plasma system in Nature beyond the measured data. That is, the predictive power of our models is admittedly and expectedly very low. Alternatively, our goal has been to demonstrate what might be achieved with a dynamical systems model if we had sufficient data to properly disambiguate the models that represent the plasma system. We did this by creating a large population of models that fit the limited data and creating and analyzing statistical distributions of their properties and their predictions. Admittedly, a large number of models giving a similar prediction do not imply that Nature is necessarily more likely to satisfy the prediction. The addition of new data, especially with alternate dilution values and frequencies of plasma measurements, could, in fact, greatly alter its predictions. However, it is, nonetheless, suggestive that a more constrained model, given superior data, might give interesting predictions that could guide further experiments. In particular, with such a more constrained model in hand, predictions could be made based on alternate dilution values and dilution frequencies as well as infusions. These predictions could be followed by laboratory experiments, and the results of these experiments could then be fed back into the model to improve it and its future predictions. Moreover, using such a constrained and improved model to guide laboratory experiments could make the search for the optimal laboratory and clinical procedure more efficient.

In order to achieve this goal, as described in Appendix A, we began by developing the structure of the DEs that would have the representational power for the plasma constituents. We considered simple decays that represented instability in the proteins in Eq. (A1) and did not require a catalyst. We also included catalyzed reactions that either created or destroyed a plasma constituent. Since catalyzed reactions were sensitive to the concentrations of each molecule attached to it, we used a model of the form Eq. (A2) for a single-species catalyzed reaction and Eq. (A4) for a reaction catalyzed by two species. With this term, the rate of production or destruction is nearly linear when the species concentrations are low compared to its threshold and nearly constant when the concentration is high compared to its threshold, and the catalyst is saturated. Additionally, we reasoned that destruction of a plasma species could only occur when the species was present in the catalyst. Therefore, we only allowed a destructive catalyst, where the coefficient is negative, when the species being destroyed was included in the catalyzed reactions. In other words, we only allowed negative coefficients in terms such as in Eqs. (A3) and (A5). This removed the possibility that the species concentrations ever became negative. All these dynamical terms were then added together to create a system of DEs of the form Eq. (A6). We also considered different classes of plasma constituents (which we call u_i , h_i , and c_i ), which might be important to the dynamical reactions creating and destroying each measured species. However, only the measured species ended up playing a role in our final population of models.

In Appendix B, we described the computational methods that we used to create a population of models that fit the data. We used a combination of a genetic algorithm, simulated annealing, and a gradient descent method. For our genetic algorithm, we began by generating a population of random models. We then trained these models using simulating annealing and gradient descent obtaining a model score that represented how close the model was to fitting the data, as seen in Eq. B2. After this, new models were iteratively created based on this population. At times these models were a combination of two models in the population, representing sexual reproduction, at other times a new model was a “mutation” of a population model with a new term or a new species in the dynamical equations, or a term or a nonmeasured species removed, and at times a brand new random model was created from scratch. With each new model, it was trained, using simulated annealing and gradient descent, and its model score was compared with that of the population. If its score was better than the worst in the population, it replaced the worst model-score model. This process was continued until the model scores of the population were satisfactorily low. Simulated annealing was accomplished by randomly choosing parameter values in a hypercube around the current point, where the size of the hypercube decreased over time. To accomplish the gradient descent method, we calculated the slope at each point and moved the parameters iteratively down the slope. The gradient descent method was done for each iteration of the simulated annealing process, and this seemed to give us the most efficient and successful results.

With the form of our DEs and our algorithms in place, we ran our program for 1,500 cpu-days to produce of population of 1,000 models as presented in the “Minimal Models” section. The model scores all fell below 0.03 as compared with the illustrative data shown in Eq. (6) and as shown in Fig.2. Using our present computational resources, we were unable to generate sufficiently large populations of models with sufficiently low scores if we increased the number of measured species much beyond this, and four species were a compromise. However, in the future, with increased access to computer clusters and with improved algorithms, it would likely be possible to scale this up to more interesting collections of data. In Fig.3, for each species, we showed distributions of the number of terms, the number of negative terms, and the number of species per term in each DE. From this initial point and throughout, we saw that there was a large difference between m ₀ and m ₂ on the one hand and m ₁, especially m ₃, on the other hand. We noted that these were the two species with the largest change to their concentrations following dilution with m ₃ standing out even further because of its great reduction. In Fig.4, we presented the fraction of single-species terms for each constituent [see Eqs. (A2) and (A3)]. In this case, we saw a strong connection between m ₀ and m ₁ on the one hand and m ₂ and m ₃ on the other hand. We pointed out that this was likely due to the fact that both m ₀ and m ₁ increased between their initial steady state to their final steady state, while the latter two decreased. The single-species reactions appear to mainly be of positive feedback loops. We displayed the coefficients and thresholds of the dominant single-species terms in Fig.5. In this case, we found that all the coefficients were positive and constructive, and that the coefficient for the dominant term in the equation for ${\dot{m}}_3$ was smaller than the others and suggested this was because it fell not only from its initial steady state but also from its diluted value, requiring less of this positive feedback.

We presented the fraction of two-species terms [see Eqs. (A4) and (A5)] for each equation in Fig.7. We noted that only five of these were present in large enough numbers in the population to make their analysis interesting. The first two were for the terms with $m_1{m}_2$ and $m_1{m}_3$ in the equation for ${\dot{m}}_1$ . Both of these terms were negative, as seen in Fig.6, and represented destructive feedback, reducing m ₁ in relation to the concentrations of m ₂ and m ₃ as well as itself. We noted that the reasons these stood out relative to other two-species terms was that m ₁ increased by the greatest amount, and, therefore, it needed these negative feedback mechanisms to stabilize it at the new higher steady state against the strong positive feedback found in the single-species term. The equation for ${\dot{m}}_3$ , on the other hand, was mostly the reverse of this. There were three terms, namely, $m_0{m}_1,m_1{m}_3$ , and $m_2{m}_3$ . The second of these was a negative feedback term, while the other two were positive (see Fig.6). It appears that the challenge of decreasing to a lower steady state required more complicated dynamics. We have already seen evidence of this in the discussion of Fig.3, where m ₃ required greater number of terms, greater number of two-species terms, and greater number of negative feedback terms. Here, we see that it needed two-species terms that are both constructive and destructive, and that achieving its dynamics is apparently a greater balance between two-species catalysts. We further pointed out that the negative feedback terms were predominantly in two-species terms, not single-species terms. We discussed that perhaps this was because negative feedback would depend on both another species concentration and its own. We pointed out that this might be why increased levels of protein concentrations appeared to be more likely in the results of the dilution experiments.1,2 Perhaps concentrations are more likely to be enhanced following dilution since reductions depend more heavily on two-species catalyzed reactions, and perhaps two-species catalysts are less efficiently explored than single-species catalysts by a stochastic evolutionary process.

In Fig.9, we exhibit the decay rates for each species and find m ₀ and m ₂ are comparable, while m ₁ is smaller and m ₃ even smaller. Once again, we suggested that this was because m ₁ and m ₃, with their larger changes following dilution, required a greater amount of their dynamics to come from balancing creative and destructive catalyzed reactions and less from breaking apart uncatalyzed.

In Appendix C, we showed a variety of examples of models from the population. We found that the population was organized into three broad groups, and we presented plots and the full DEs for four models from each group. The first group was composed of models, where $m_0,m_1,\text{and}{m}_2$ began increasing immediately following dilution while m ₃ immediately began decreasing. Plots for exemplary models were shown in Fig.15, while their DEs were shown in Eqs. (C1) through (C4). The second group contained models where all of $m_0,m_1,m_2$ , and m ₃ began increasing immediately following dilution. We presented plots for illustrative models in Fig.16 and their accompanying DEs in Eqs. (C5) through (C8). The last group was the same as the first group except that the concentrations oscillated up and down before approaching their final steady states. Example models were presented in Fig.17, and their DEs were given in Eqs. (C9) through (C12).

FIG. 15.

A sampling of model dilution plots from the “Minimal Models” section from Class 1 (see the beginning of this appendix). The system of differential equations for these plots can be found in (a) Eq. (C1), (b) Eq. (C2), (c) Eq. (C3), and (d) Eq. (C4).

FIG. 16.

A sampling of model dilution plots from the “Minimal Models” section. The system of differential equations for these plots can be found in (a) Eq. (C5), (b) Eq. (C6), (c) Eq. (C7), and (d) Eq. (C8).

FIG. 17.

A sampling of model dilution plots from the “Minimal Models” section. The system of differential equations for these plots can be found in (a) Eq. (C9), (b) Eq. (C10), (c) Eq. (C11), and (d) Eq. (C12).

In the “Dilution Variations” section, we explored and illustrated the properties we might learn about the dynamical system and about alternate dilution values if we had greater data and a more realistic, more robust model. In Fig.10, we considered the behavior if we diluted all constituents by a flat amount but at different dilution values. We found that for any dilution values, the system was likely to end up in the OFS, but if, instead of a dilution, an infusion of all four species were achieved, the system would evolve to an AFS. Although it is not likely we would infuse all the constituents of plasma in this manner, the AFS in this plot showed up in every other case we discussed in this section. In Fig.11, we considered what would happen if we diluted all species at the original flat 50% except for one species, which was diluted at a variety of levels. If the single species diluted at a different level were one of m ₀, m ₁, or m ₂, then the system would typically still evolve to the OFS, even if the single species was instead enhanced. However, the model proved especially sensitive to the dilution value of m ₃. Indeed, we found that any dilution of m ₃ above approximately 70% would lead the system to finish at the AFS. Interestingly, and as mentioned, this AFS is the same as in Fig.10. We showed a similar behavior in Fig.12, when all the species were diluted at the original flat 50% dilution except for two species, which were diluted independently at a range of values. The region was again split into one region where the system evolved toward the OFS, while in the other region, the system culminated in the same AFS as before. Once again, we saw that the system was most sensitive to the dilution of m ₃. Just as in the “Minimal Models” section, we see that the model is much more sensitive to m ₃, presumably because it decreases by the greatest amount from its initial steady state to the OFS, which leads to greater sensitivity in its DE.

We also considered what would happen if we held all the species at a fixed 100% value while diluting only one of the species at a variety of levels and showed the results in Fig.13. Interestingly, we found that if we diluted or enhance any single species in the direction of its OFS, that is where the system evolved, while if we diluted any species in the opposite direction of the OFS, it would instead evolve to the same AFS as before. We further noted that this study showed that the initial steady state was unstable to small perturbations, and that we considered this to be an obvious shortcoming of this model. However, we discussed that artificially implementing a region of stability around the initial steady state would be arbitrary, and that, instead, we should seek experimental dilution data within the stability region as well as outside in order to guide our models. In Fig.14, we considered holding all but two species fixed at the initial steady state value and independently varying the dilution of the other two. We found a similar behavior to Fig.13, where diluting in the direction of the OFS leads the system there and diluting in the other direction leadsto the AFS. We also noted that the presence of an AFS and the near universality of it among these models, even though it was neither part of the training data nor explicitly programed, suggest that, with better data, we might predict ASSs which might be useful in the lab or the clinic and dilution regimens, which might achieve them. In Appendix D, we presented various dilution scenarios for a single model for illustration.

As we have stated repeatedly throughout this note, the available data are not yet sufficient to make robust predictions, and, therefore, the specific predictions made in this note should be taken as nothing more than illustrative of what sort of predictions we might make if we gathered better data. In particular, the data that would be, the most helpful would be:

Hourly measurements before dilution to determine the circadian cycle for the plasma constituents. Thebehavior post-dilution could be impacted bythe timing of the dilution events since they might be higher or lower in the morning relative tothe evening. This would also give a more realistic model and understanding of the plasma system.

Hourly measurements following dilution. Such datawould have a large impact on disambiguating the various models that are only fit to the endpoints since it would constrain the models to fit the data in between the initial and final steady states.

Measurements with a variety of dilution amounts, including some infusions. This would further disambiguate the models, establish a region of stability around the steady states, and lead to greater predictive power for new experiments.

We encourage the experimental community to carry out these experiments. We would also like to encourage the release of more precise concentration levels for each constituent, which would allow analysis beyond merely stating whether a concentration is relatively high or low. This could be important, in particular, in better understanding the catalyzed dynamics and their threshold values.

Authorship Contribution Statement: DC: Visualization and Writing. NC: Conceptualization, Methodology, Software, Formal Analysis, Investigation, Writing, Visualization, Supervision, and Project Administration.