Modeling the Effects of Multiple Myeloma on Kidney Function

Multiple myeloma (MM), a plasma cell cancer, is associated with many health challenges, including damage to the kidney by tubulointerstitial fibrosis. We develop a mathematical model which captures the qualitative behavior of the cell and protein populations involved. Specifically, we model the interaction between cells in the proximal tubule of the kidney, free light chains, renal fibroblasts, and myeloma cells. We analyze the model for steady-state solutions to find a mathematically and biologically relevant stable steady-state solution. This foundational model provides a representation of dynamics between key populations in tubulointerstitial fibrosis that demonstrates how these populations interact to affect patient prognosis in patients with MM and renal impairment.

Thus far, mathematical modeling related to MM has focused on the breakdown in bone remodeling process caused by malignant plasma cells [4][5][6] . In this paper, we focus instead on the kidney damage that occurs in some patients with MM because of the significant effects of kidney dysfunction on patient prognosis.
Several studies have reported inferior overall median survival time for patients who present with renal impairment, which occurs in approximately 50% of patients with MM 7 . A study by Knudsen et al. in 2000 grouped patients according to degree of renal impairment: none (defined as creatinine <1.47 mg/dL), moderate (creatinine 1.47-2.26 mg/dL), and severe (creatinine >2.26 mg/dL). They found a median survival of 36 months for patients in the first group, 18 months for patients in the second group, and 13 months for patients in the third group. Furthermore, the study concluded that reversibility of renal failure slightly improved prognosis 8 . More recently, a 2014 retrospective study that included patients from 15 Swedish hospitals used estimated glomerular filtration rate (eGFR) to study effects of renal impairment on overall survival time with two types of treatment. Overall median survival time (33 months) for patients presenting with renal impairment (defined as eGFR < 60 mL/min/1.73 m 2 ) was significantly less than overall median survival time (52 months) for patients with no renal impairment 9 . A 2015 study by Gonsalves et al. found that the median survival rate for patients with renal insufficiency (creatinine >2 mg/dL) was 42 months while median survival time for patients with normal renal functions was over twice that, at 99 months 10 . Patients in this study were treated with one or more novel agents, and it was found that patients who demonstrated improved renal function after treatment had improved prognoses compared to patients with no response. However, the improved survival rates were significantly lower than the prognosis of patients with normal renal function. This disparity between expected survival time in patients with and without renal insufficiency is the motivating clinical problem for our model.
We present a continuous and deterministic model of dynamics between the proximal tubule cell population and free light chains in the kidney of a healthy patient and then a model that incorporates the effect of malignant plasma cells in patients with MM. We analyze our model for steady-state solutions to find a mathematically and biologically significant stable steady-state solution and find computational results for the model incorporating malignant plasma cells. Our model captures the main players and relationships in the process of renal interstitial fibrosis caused by MM and is meant to translate some basic understanding of MM for future work on predictive and prognostic models useful for patient-specific medicine.

Renal Function
In order to examine the effect of MM on the kidney, it is necessary to first consider normal renal function. Humans have two kidneys, made up of systems of tubules, called nephrons, that are the working units of the kidney. Each kidney has approximately 1 million nephrons (Fig. 1). The kidneys are the body's filtering system, and help maintain the body's homeostasis. The kidneys keep substances in the human body in balance by regulation and removal of metabolic waste. The main functions of the kidney include regulation of water and electrolyte balance, which affects arterial blood pressure. The kidneys also monitor the excretion of hormones and the blood levels of prescription drugs that affect body function, as well as regulating red blood cell production and vitamin D production. Most of the kidneys' work involves transportation of water and solutes between the blood and filtration systems in the kidney. Any substance not transported back to the blood is excreted in urine.
Renal disease in patients with MM is usually present as proteinuria (excess serum proteins in the urine) 11,12 . In general, kidney diseases related to MM result from the kidneys' reduced ability to properly filter substances. There are two types of common tubular kidney damage that we will consider: proximal tubular cell injury and cast formation. The sources of this damage are malignant plasma cells, which produce monoclonal proteins called monoclonal immunoglobulin (Ig) light chains (also called M proteins). Ig light chains are subunits of antibodies, which are produced by plasma cells to fight infection, and antibodies express one of two types of light chain: kappa light chain or lambda light chain. Monoclonal myeloma plasma cells will overproduce one specific type of Ig light chains. The most common cause of severe renal failure in patients with MM is tubulointerstitial pathology resulting from high circulating concentrations of monoclonal Ig light chains 11 .

Proximal Tubule Cell Injury
Proximal tubule cell (PTC) injury is caused by free light chain interaction with PTCs. Free light chains can be toxic to PTCs by blocking transport of glucose, amino acids, or phosphates, activating tubulointerstitial fibrosis, and causing excess free light chain endocytosis 13 . Endocytosis is the process by which cells absorb proteins by engulfing them. It is important to note that not all monoclonal free light chains are toxic to the kidneys, and in fact, many patients with high amounts of serum free light chain have no renal impairment 7 . It appears that toxicity depends on the structure of a particular individual's free light chains' 3D structure or protein folding 11,14 . While high amounts of light chains can be a sign of MM, a more useful laboratory parameter is the serum kappa to lambda ratio. Normal kappa to lambda ratio is 0.26-1.65, compared to 0.37-3.1 in situations of renal impairment 15 . When the level of either kappa or lambda light chains is very high and the other is normal to very low, the ratio is considered abnormal, which suggests the presence of clonal plasma cells.
One way free light chains can be toxic to PTCs is by activating tubulointerstitial fibrosis. Tubulointerstitial fibrosis is the process initiated by the interaction between proximal tubule cells and free light chains, which activates inflammatory pathways in the kidney. Sustained inflammation causes the excess accumulation of extracellular matrix (ECM), which is eventually replaced by scar tissue 16 . ECM, which is made up of proteins and collagens, provides structural support to surrounding cells and most cells cannot survive unless they are anchored to the ECM. The scar tissue that replaces ECM is part of the formation of excess fibrous tissue that characterizes fibrosis. This process is considered to be largely irreversible, and leads to the loss of function of proximal tubule cells and end stage renal disease (ESRD).
Because tubulointerstitial fibrosis begins with the interaction between proximal tubule cells and free light chains, the main goal of treatment is to reduce light chain production by killing the malignant plasma cells 12,17 . Initial treatment includes chemotherapy drugs, hydration, and in some cases the use of bisphosphonates to lower calcium levels. Also used as treatment is plasmapheresis, which involves removing the blood plasma from the body, treating it (reducing plasma concentration of light chains), and returning it to the body. This is similar to dialysis, which is used to remove waste from the blood. Both plasmapheresis and dialysis use machines to perform the kidneys' usual job of filtering the blood. Renal transplantation is usually not considered because of the poor prognosis of patients with MM 12 .

Flowchart
To better identify the various cell populations involved in the processes we model, we present a flowchart in Fig. 2.
The increased monoclonal protein production by the myeloma cells leads to an increased level of free light chain molecules circulating in the blood. These free light chains are either endocytosed, or precipitated. In the primary situation that our model addresses, the increased free light chain production leads to increased light chain endocytosis by proximal tubule cells via cubilin/megalin complex 12 . Cubilin and megalin are two endocytic receptors that play important roles in renal tubular clearance and reabsorption of proteins. They initiate receptor-mediated endocytosis, a process by which cells internalize molecules. This involves an inward budding of plasma membrane vesicles containing the monoclonal proteins with receptor sites specific to the molecules being internalized. This increased light chain endocytosis activates NF-κ B and MAPk in the proximal tubule cells. NF-κ B is a protein complex involved in regulating the immune system's response to inflammation, and is responsible for cytokine production. Mitogen-activated protein kinases (MAPk) direct the cellular response to mitogens and proinflammatory cytokines.
The activation of NF-κ B and MAPk initiates the production of several different types of cytokines and growth factors by the proximal tubule cells: IL-6, CCL2, IL-8 and TGF-β. IL-6 is secreted by T-cells and macrophages to stimulate immune response, and acts as a proinflammatory cytokine. IL-8 is produced by macrophages and epithelial cells. CCL2 recruits memory T-cells and dendritic cells to inflammation sites. TGF-β is a protein that controls cell growth, apoptosis and proliferation. These cytokines and growth factors initiate proinflammatory and fibrotic pathways, and initiate Epithelial-Mesenchymal Transition (EMT), type 2. During EMT type 2, polarized epithelial cells (such as those that line the kidney tubules, in our case, proximal tubule cells) change to assume mesenchymal cell characteristics. This allows these cells increased migratory ability to migrate to an infection site, increased resistance to apoptosis, and increased production of ECM material. This all plays a part in renal interstitial fibrosis, the sustained inflammation in proximal tubule epithelial cells. Fibrosis causes a disruption in the normal genesis and breakdown cycle of ECM, which leads to excessive ECM accumulation 18 . Eventually, scar tissue replaces ECM accumulation, and causes loss of function of PTCs. Ultimately, end-stage renal failure can develop.
In the secondary situation in our flowchart, non-endocytosed free light chains precipitate, forming solids called tubular casts within the kidney tubules. These casts are formed by the reaction of Ig light chains with Tamm-Horsfall protein. The casts partially or totally block the kidney tubules, which increases intraluminal pressure, reduces glomerular filtration rate (GFR), blood flow, and tubular clearance of the light chains, which increases serum light chain levels, resulting in a never-ending cycle. Unless the casts are removed, the result is permanent nephron loss.
Current kidney physiology modeling focuses on modeling chemical exchange between compartments in the kidney, and on modeling GFR [19][20][21] . GFR depends on age, sex, and body size, and gives a good indication of how well the kidney is functioning and filtering substances in the body. To our knowledge, there is no known prior mathematical work in modeling the above process of renal tubulointerstitial fibrosis caused by MM.

Model Development
To create our mathematical model, we use modified power law approximations, developed by Savageau and Voit 22,23 . Power laws are useful here because they capture the nonlinearity specific to biological systems such as this one, but are comparatively easy to work with analytically. Power laws have the following form: where the first term represents growth of the X i population affected by X j populations, and the second term represents death or removal of the X i population affected by X j populations. The parameters α k are growth or proliferation rates and the parameters β k are death or clearance rates. Based on the biological background from Fig. 2, we focus on the populations of PTCs, FLCs, and renal fibroblasts for our initial model for normal dynamics, and then include the tumor cell equation for our model that simulates dynamics in a patient with MM and renal degradation.
In our simplified model of normal dynamics, the growth of PTCs is governed by its own proliferation rate and the population of PTCs decreases only through apoptosis. The growth of FLCs increases at a natural production rate and decreases by a natural renal clearance rate. The growth of renal fibroblasts increases at a natural production rate and decreases by apoptosis.

Model of PTC and FLC Dynamics in the Kidney without Tumor.
Using the biological background and power laws discussed above, we construct a system of ordinary differential equations (ODEs) for the PTCs and FLCs in the kidney of a healthy patient: where we define (x) + = max(x, 0) and P(t), L(t), and F(t) are the populations of PTCs, FLCs, and renal fibroblasts. We use the maximum function to preclude decrease in any of the populations in the situations L > L S , >L L, and >F F. Because these are growth terms, it does not make sense to have a negative value (the decay terms capture decrease in these populations). The use of the maximum function controls the growth to allow a minimum of 0, or no effect, and maximum of 1. In this system, Term A represents the proliferation rate of PTCs, where the proliferation rate of PTCs decreases in the presence of increased FLCs. Note that in a healthy patient the level of FLCs will stay approximately constant. The parameters β P is the PTC proliferation rate, L S is the saturation level of FLCs, and g 1 and g 2 describe the strength of FLCs on PTC growth and the strength of PTCs on its own population growth, respectively. The use of g 1 and g 2 in this equation reflect the use of modified s-form equations and that the growth of the population of PTCs is not dependent on a true logistic equation. Term B represents a decrease in the population of PTCs due to apoptosis. In the FLC equation (2), Term C represents the increase in number of circulating FLCs, where we define , L min is the minimum number of circulating FLCs and γ L is the growth rate of FLCs. This entire term is defined to be positive or zero. Term D represents the natural clearance or removal of FLCs where μ L is the natural clearance rate. In the renal fibroblast equation (3), Term E represents the increase in number of fibroblasts where we define and F max is the maximum percentage of renal fibroblasts. As with Term A and Term C, this term is defined to be non-negative. Term F represents loss of renal fibroblast population due to apoptosis. In the FLC equation (2) and the fibroblast equation (3), we assume linearity in the growth of Term C and Term E, and therefore the exponents on these terms are set to 1 for simplicity.

Model of PTC and FLC Dynamics in the Kidney with Tumor. To simulate the dynamics between cells in the kidney in a patient with MM and renal degradation, we add three terms to our previous system, as well as incorporate a fourth equation representing the growth of malignant plasma cells via tumor equation T(t).
We chose to use the differential equation form of the Gompertz model to model tumor density. The Gompertz model assumes a time-dependent growth rate and a carrying capacity for tumor load, and results in a sigmoidal shape 24,25 . It has been previously used to simulate the growth of malignant plasma cells in patients with MM 4 . The new system of ODEs is shown below:   4 . All other parameters were obtained for the model using heuristic parameter estimation to a generic response. L min = 10 mg/L represents the expected amount of FLCs in any patient, and L S = 1000 mg/L is the threshold at which renal impairment is typically observed.
Scientific RepoRts | (2019) 9:1726 | https://doi.org/10.1038/s41598-018-38129-7 T S P As in the previous model, we define In this set of equations, Term I now represents the loss of PTC due to EMT (which is precipitated by the excess FLC population). The parameter γ F is the fibroblast growth constant, demonstrating that as the number of fibroblasts increases, more PTCs have undergone EMT. The exponents g 3 and g 4 represent the strength of FLCs on PTC transition and the strength of PTCs on their own transition via EMT. In (5), the new term is Term K, which represents the increase in the number of FLCs produced by malignant plasma cells. The parameter γ L represents the growth rate of these plasma cells. The exponents g 5 and g 6 represent the strength of the tumor on FLC growth and the strength of FLC on its own growth. In (6), Term N is the same term as Term I, but is a source of growth in the fibroblast population. The final equation (7) is the Gompertz tumor growth equation, where T S represents tumor carrying capacity, or maximum spread of the tumor.
The parameters for the Gompertz model (γ T , T S ) are scaled to a dimensionless model, e.g. maximum tumor load is set to be 100. The tumor growth parameter values were taken from published work by Ayati et al. 4 . All other parameters were obtained for the model using heuristic parameter estimation to a generic response, provided as a human-based aggregate of observation by the clinical author. Because of the invasiveness of kidney biopsies, obtaining direct information about quantities of renal PTCs and fibroblasts in particular is difficult, and little data is available for rigorous parameter estimation techniques. Validation of this model relies on further data. The model is a foundation for incorporating clinical data into a quantitative predictive system, which remains future work. The parameters used to obtain the initial, illustrative results in this paper are listed in Table 1.
(P * , L * , F * ) Requirements for Stability Table 2. Summary of the results of steady-state analysis on the model described by Equations (1)-(3). Because we expect a healthy patient to have non-zero amounts of PTCs (P), FLCs (L), and renal fibroblasts (F), only one of these eight steady-state solutions is both mathematically and biologically significant:  Table 2 with their requirements for stability based on analysis of the Jacobi matrix. We expect that healthy patients have non-zero levels of all three populations of interest, so only one of these steady-states is both mathematically and biologically significant: The requirements for stability of this steady-state solution are γ F > μ F , γ L > μ L , g 2 < 1. This indicates that the growth rate of the fibroblast population, γ F , must be greater than the fibroblast death rate, μ F . Similarly, the growth rate for the free light chain population, γ L , must exceed the clearance rate μ L . Finally, g 2 < 1 indicates that the strength of PTC on its own growth must be nonlinear (more specifically, less than linear). To estimate initial parameters for the model, we used known values where available and scaled the rest to maintain a stable steady-state when P * = 100 for the case above; that is, a healthy patient who starts with 100% of PTC population will remain healthy. In a healthy patient, normal FLC levels are approximately 10 mg/L, so we set L min = 10 and L S = 1000, as kidney problems typically occur at FLC levels above 1000 mg/L 26 . We also analyzed the stability of the steady-state (P * , L * , F * ) = (100, 10, 100) computationally using the ode15s solver in MATLAB. Results of this analysis are included in Fig. 3, demonstrating that for (P * , L * , F * ) = (100, 10, 100) our model has a stable steady-state solution and in a patient with no myeloma tumor cells present, any perturbation in normal conditions will result in eventual return to homeostasis conditions (and the steady-state solution) as expected.
Computational Results for Model with Tumor. Figure 4 shows computational results generated using MATLAB ode15s, with initial conditions P(0) = 100, L(0) = 10, F(0) = 100 and T(0) = 1. Note that these initial conditions represent that the patient has 100% of PTCs, normal levels of FLCs at 10 mg/dL and 100% of fibroblast  7)) using the parameter values in Table 1, for initial conditions P(0) = 100, L(0) = 10, F(0) = 100 and T(0) = 1, over a time period of 160 days using the model given by Equations (4)- (7). This simulation represents a fast-acting tumor that is not treated. We see that as the percentage of tumor cells increase in (d), the amount of free light chains increases beyond the threshold amount of 1000 mg/L in (b) and amount of renal fibroblasts begins to slowly increase in (c). Decline in proximal tubule cells as a response to (b-d) is demonstrated in (a). population. These three values all are consistent with a healthy patient with no myeloma cells present in the body. However, the initial tumor density begins at 1% and γ T = 0.05 days −1 , representing a quickly growing population of malignant plasma cells. These results are consistent with the biology described in Fig. 2: as the tumor cell population grows, the free light chain population increases, the proximal tubule cell population begins to decrease and the fibroblast population increases. Figure 5 demonstrates the relationships between two different sets of variables by plotting the PTC population against the FLC population in 5 (a) and the FLC population against the myeloma tumor (or monoclonal myeloma protein) population in 5 (b). Figure 5(a) shows results for the model run over 160 days. As the FLC population increases, the PTC population decreases. When the FLC population reaches 500 mg/L, a common threshold for patients with kidney damage, the PTC population has decreased to between 40-50% of the original population. Figure 5(b) demonstrates that as the tumor population increases, there is a delay in the response of increase to the free light chain population.
The computational results for this model simulate the expected biological relationships between the key populations involved in renal impairment in patients with MM and lay a foundation for a model incorporating the effects of treatment. Because of risks involved, patients presenting with severe renal impairment are often excluded from treatment trials 7 . An in silico experiment via computer simulations of a mathematical model could give researchers more opportunities to investigate the impact of renal impairment on overall median survival time and patient prognosis.

Conclusions
We have developed a mathematical model for the major cell and protein populations involved in proximal tubule cell injury due to free light chains produced by myeloma cancer cells. The system of ODEs captures the qualitative behavior of the cell populations based on a biological understanding of the process taking place inside the kidneys, that is, that the proximal tubule cells in the kidney undergo apoptosis as a result of increased production of free light chains by myeloma cells and inflammatory response involving renal fibroblasts.
This model is a mathematical representation of the processes that take place in patients with MM that present with renal injury at diagnosis. Several simplifications were made when developing the model. These include the decision to focus on only four populations within the body and hence on parameters that reflect the effects of those populations on each other, but not other influences. In addition, this model does not take into account other factors that may influence renal function including other co-existing diseases such as diabetes mellitus and hypertension or medications.
Further work in this area includes parameterization of the theoretical model utilizing relevant markers of renal function and myeloma burden from patients suffering from this disease. It is our hope that this model will lead to further models, incorporating essential processes as needed, that could be calibrated with training sets from clinical data and then validated by comparing model predictions against actual patient outcomes. Increasingly faithful models could be used to further investigate the relationship between renal function and serum FLC levels at presentation, or to predict the likelihood of renal function recovery following myeloma therapy for those patients presenting in renal failure. Even for patients without overt signs of renal failure it has been assumed that upon diagnosis of myeloma, some degree of renal injury has already been sustained. These modeling efforts could  7)) using the parameter values in Table 1, for initial conditions P(0) = 100, L(0) = 10, F(0) = 100 and T(0) = 1, over a time period of 160 days. The proximal tubule cell and tumor populations are interpreted as percentages, and free light chain population units are mg/L. This simulation represents a fast-acting tumor that is not treated, and shows the relationship between proximal tubule cells and free light chains, and the relationship between tumor growth and free light chain growth, respectively. lead to a better understanding of this sub-clinical damage and therefore have an impact on clinical strategies which could be employed to ensure optimal renal function. Success would provide clinicians with a valuable tool with genuine prognostic and predictive capabilities.