Mathematical and numerical modelling of copper transport in yeast

The transport and regulation of metals in eukaryotic cells is a complex process, dependent on protein transporters that respond to cell needs. The application of dynamic mathematical models can provide valuable insights into these transport mechanisms. Mathematical simulations of transport processes may not directly predict transport mechanisms but can guide experimental design or identify inconsistencies between observation and hypotheses. Copper is an essential metal in eukaryotic cells as a catalytic co-factor in metallochaperone proteins and is therefore tightly regulated in living systems, making it valuable for quantifying biological transport mechanisms. In order to test our modeling system, a culture of baker’s yeast (Saccharomyces cerevisiae) was grown, copper concentrations were obtained from the cells and growth media, and a mathematical model was developed to investigate transport mechanisms between the growth media and the cells. A model based on conservation of mass was presented as a system of equations upon which to develop. This system of equations was developed to include an active transport term that describes a homeostatic concentration that cells actively maintain through negative feedback, and with a delayed activation, the model was more accurate at predicting the experimental data. The hypothesis and dynamic model derived in this work provide a novel framework that may be applied to additional metals or used to describe other transport mechanisms in biological systems.


Introduction
Copper is an essential element to support metabolic processes in the body [1]. Existing as Cu(I) or Cu (II), copper participates in several biologically-mediated oxidation-reduction reactions where the atom transports electrons. It is critical that the amount of copper is carefully regulated, as the presence of excessive copper can result in damage to delicate intracellular organelles by oxidative stress [2]. To prevent this risk, sophisticated strategies for maintenance of copper homeostasis were selected for during evolution of eukaryotes. For example, human metallochaperones, like Menkes (ATP7A) and Wilson (ATP7B) proteins, safely bind and transport metallic copper through the cell [2]. Copper was selected as the model metal in this work due to its role in a number of ailments such as Menkes' and Wilson's disease, which occur due to mutations in the genes coding for both of these metallochaperones, ATP7A and ATP7B, respectively [3]. It has been demonstrated that ATP7B affects the stable isotope composition of copper in blood serum [4] and antibiotic treatment can further alter the isotope composition in gut flora [5]. Moreover, copper is an important element in biomedical applications such as copper-free click chemistry [6] and copper-based antitumor complexes [7].
Because of the vital nature of copper in living systems, the investigation of processes involving this element is an active area of research. One opportunity to expand our knowledge of copper interactions is the development of computer models to model dynamic copper exchanges such that changes in the distribution of copper can be visualized using the principles of inverse modelling. There are many contemporary uses of inverse modelling, such as computed tomography, where data and models are integrated to create images of directly unobservable areas. For instance in the modelling of tides [8], seismic modelling [9], the interpretation of airborne electromagnetic data [10], or magnetotelluric analysis [11]. The challenge with this approach is that real systems are dynamic and the processes responsible for copper transport and exchange are continuously evolving and changing in time. Therefore, a simple model based on a single mechanism is unlikely to produce results that reflect 'real-world' conditions. Further, the actual processes may not be known, and it is necessary to create empirical models that predict the behaviour without detailed knowledge of the underlying chemical reactions. In this paper, we report on the development of a computational model to describe the movement of copper from a nutritional solution into a growing yeast culture. Beginning from the simple assumption that copper incorporation is driven by diffusion into the cells, the model was expanded to include more dynamic terms that take into account the actual number of cells and the stage of development of the culture. Experimental data obtained from sampling a culture of baker's yeast, Saccharomyces cerevisiae, were the basis for the model and the quality of the fit was assessed using numerical root-finding algorithms; a process known as inversion. Baker's yeast is a simple organism that well models the behaviour of human cells.
This paper describes the evolution of the model from a simple diffusion mechanism, to the inclusion of terms that describe an active transport process, and finally the addition of a non-linear sigmoid transition function, which was necessary to simulate the lag in the uptake of copper as the culture acclimated to its growth media. The final form of the model is a framework upon which the transport of copper (and other biologically important metals) can be modelled and investigated.

Modelling and results
Population modelling Yeast was grown in the presence of glycerol as a non-fermentable carbon source to maintain the culture respiring throughout all phases of growth. During the growing period, the culture was sampled at regular intervals. The samples of the culture were centrifuged to isolate cells from the growth media. The resulting media and cells were analyzed for copper concentration.
The growing yeast culture was modelled using its growth rate r, which was defined as the change in the number of cells Y(t) per unit time. The carrying capacity K was defined as the maximum number of cells in the culture. The logistic model was frequently used as a model for population growth when there is a carrying capacity and the population is not capable of growing exponentially. The logistic model is essential for modelling the growth of yeast in a finite growth medium and forms the basis of this work. The logistic model is represented by the nonlinear differential equation, such that when Y=K, the rate of change of the number of cells Y(t) became zero and the culture ceased to increase in number. The solution to equation (1) was the analytical form of the nonlinear sigmoid curve, where Y 0 is the initial size of the yeast culture. The cell density Y(t) was estimated as the total number of cells obtained from measurements of the optical absorbance of the culture. The analytical solution, equation (2), was fit to the number of cells per unit volume sampled along the growth curve to obtain estimates of the growth rate parameter r and carrying capacity K (see figure 1). Fitting the number of cells in a unit volume through nonlinear regression to the analytical solution of the logistic equation returns a growth rate parameter r of 0.41 h −1 with a standard deviation of 0.01 h −1 and a carrying capacity K of 5.02×10 8 cells ml −1 with a standard deviation of 0.06×10 8 cells ml −1 .

Mass balance modelling
With knowledge of the growth rate parameter and carrying capacity, it was possible to focus on modelling the copper exchange between the growth media and the cells in the growing yeast colony. The model describes a small change in the concentration of copper in a cell ∂N c over a small increment of time, dt, Where (∂N c,I -∂N c,O ) represents the net uptake of copper by cells through the membrane and ∂N c,Y accounts for redistribution of copper during population growth. A first-order Taylor series expansion was used to approximate the post-mitosis cellular concentration as the sum of a pre-mitosis concentration and ∂N c,Y, due to the increased number of cells,

( )
Considering the copper mass redistribution during population growth first, mass must be conserved during mitosis, so, Note that during population growth Y'/Y>1 so N' c <N c after cell division, as one would expect. The equation was rearranged to isolate for the post-mitosis concentration of the cell N' c , From substitution of equation (6) with equation (4), we have,

( )
Using another first order Taylor expansion for the post-mitosis population, this can be written as It is important to note that ∂N c,Y was negative when the population was growing because copper was redistributed among each subsequent generation yielding a smaller average concentration per unit biomass. Now consider the influence of copper uptake through the cell membrane (∂N c,I -∂N c,O ) for which several models were possible. The simplest model for copper uptake was based on diffusion, Where the first term described changes to copper concentration per cell driven by a concentration gradient between growth media and cytosol. It was clear that a diffusion model could not fit these data since the maximum value of N c =N g and a diffusion model predicts that dN c /dt≈0 when N g ≈N c , and the data did not support this. In fact, the uptake of copper by cells was negligible compared to the available copper in the growth media N g , so it was reasonable to assume, Equation (11) describe a transport process where copper diffuses from the growth media into the cell. Inversion of the data with the model developed in equations (11) provided no meaningful solution because diffusion was not the dominant method utilized by cells to transport copper. A different term must be included in the model to establish an amount of copper that the cells actively regulate. This active transport of copper from the growth media to the cell occurred independently of concentration gradients. The incorporation of the active transport mechanism in the model equation required an additional term to fully describe the movement of copper to the cell. An active transport model for cellular copper uptake was proposed, in which the rate of uptake was proportional to the concentration in the growth media N g with rate parameter r a and the cellular copper concentration N c relative to a steady-state value N c,K regulated through a negative feedback mechanism. We determined the best fit values of r a and N c,K through nonlinear regression (see Materials and Methods section for details), Figure 3 displays the results of the inversion for the active transport process. The rate parameter r a was estimated as 5.0×10 -5 h −1 and the steady-state regulated term N c,K was 0.28 fg of copper per cell. Perturbations of 20% were applied to parameters r a and N c,K to observe how the stability of the model varied with respect to uncertainty in each parameter over a range of values. The model was more sensitive to changes in N c,K than to changes in r a and this is understandable. The parameter r a controls the rate which the culture reaches its target concentration. A rate greater than 5.0×10 -5 h −1 will result in a fast accumulation of copper, and a lower rate will result in a slower accumulation. Regardless, this model is deficient because it is defined by its unrealistically rapid transport of copper. Lowering the rate parameter r a to match the delayed response of the culture will require an increase in the target concentration N c,K to fit the data. This is not possible because this concentration term was determined experimentally and it follows that this model is erroneous.
The model, based on a steady-state regulated copper concentration, improved on a diffusion model by accurately capturing the maximum concentration N c , but failed to account for the lag in copper uptake at the beginning of the population growth. This lag time, attributed to the period of time cells adapt to the fresh medium, could be accounted for by introducing a sigmoid function 1 / (1+e -rs(t-t0) ) with two additional parameters r s and t 0 which accounted for the activation rate and the lag time before copper uptake activation. Again, the best fit values for r a , N c,K , r s , and t 0 were determined through nonlinear regression.
The uptake rate parameter r a was estimated with a rate of 3.0×10 -5 h −1 . The steady-state regulated term N c,K was estimated again as 0.28 fg per cell and delay parameter t 0 was set to 9.5 h. The activation parameter r s was estimated as 1.9 h −1 . Perturbations of the parameters were performed to assess the sensitivity of each parameter in the model. Varying the parameter t 0 adjusts the offset as predicted so that a lower value for this parameter meant that copper uptake occurred earlier but it did not affect the slope. Perturbations of 20% on the activation parameter r s affected how abrupt the activation was, with lower values resulting in a more gradual uptake of copper and higher values (representing an activation over a shorter period) yielding sharper transitions. Of the three models presented, this system of equations (14), and (1) yielded the best fit with the data through inversion (see Materials and Methods section for more information on the algorithm used).

Discussion
Though a basic diffusion model could simulate the absorption of copper into cells, there was no method of limiting this uptake, and diffusion would continue until a concentration balance was achieved between the growth media and the cells such that dN c /dt≈0 when N g ≈N c . This process did not consider that cells actively regulate intracellular copper concentrations through a negative feedback mechanism [12].
Therefore, the diffusion model was rejected and replaced with a set of equations that included a term that described copper transport by enabling cells to absorb a certain amount of copper for biological functions to proceed effectively and then actively participating in regulating copper to maintain this level. The model including this regulated term N c,K was developed from the hypothesis that cells would achieve a stable steadystate copper concentration that was mediated through biological action based on the needs of cells as opposed to passive diffusion into the cell. The inclusion of the active process improved the fit of the model to the experimental data. However, the system of equations was still inadequate in modelling the lag in copper uptake observed until around t=9 h of the experiment. Incorporating a sigmoid function with a delay parameter t 0 and coefficient r s modelled a gradual activation of the transport process. This sigmoid function emulated the lag phase experienced by yeast as it adapted to its new growth media. Upon activation of the transport function, the uptake of copper was rapid and quickly reached a steady-state equilibrium and remained at equilibrium as cells uptake or effused copper, but this transition was more gradual with a sigmoid function.
After fitting the model to the data, it was clear that sampling the culture over the full experimental time frame would be important in future experiments. The number of cells was low at the start of the experiment and not much copper was uptaken by the yeast biomass. Furthermore, the model developed in this study focussed on the bulk transport of copper between the growth media and cells but did not consider the intracellular transport of copper-containing metallochaperone proteins. It would be interesting to examine the intracellular distribution of copper among different organelles. For example, a non-proteinaceous pool of labile copper in the mitochondrial matrix of yeast is known to be essential for assembly of cytochrome c oxidase [13,14]. Therefore, subcellular fractionation and purification of mitochondria could aid in investigations aimed to monitor the intracellular transport between the labile copper pool predicted by Fahrni [15] in the cytosol and this key organelle, where copper concentrations have been estimated to reach 200 μM under respiratory growth conditions [14,16]. The yeast S. cerevisiae has been widely used to investigate copper metabolism [17]. A collection of mutants lacking non-essential copper transporters, cuproproteins or metallochaperones is available. Combining our sensitive copper detection method and model developed herein with these mutants represents a unique powerful high throughput approach to dissect how cells maintain copper homeostasis and how deviations result in disease. For example, bioavailable copper is mainly incorporated into yeast cells via the high-affinity transporter Ctr1, but deletion of Ctr1 does not appear to affect copper uptake, although cells are much more sensitive to excess copper [18]. Comparing copper concentrations in wild type versus cells lacking Ctr1 over time could provide profound insight into the effects of unregulated copper uptake and toxic copper thresholds. The investigation of metallochaperone transport at a subcellular level is an important continuation of this work because this information could provide a deeper understanding of how metals are transported within the cell and the role metals play in disease.
This work utilized IDMS to obtain copper concentration data through mass spectrometry but there is the opportunity to investigate the isotope abundance fractionation of the two stable isotopes of copper, 63 Cu and 65 Cu. Equilibrium-based isotope fractionation was predicted due to energy differences in metallochaperone copper binding sites [4,19] suggesting that there are distinct isotope abundance variations between principal organelles involved in intracellular copper regulation.
The model detailed in this work provides insights into the process of copper uptake while avoiding detailed knowledge of the actual mechanism involved. This model functions as a framework for this research to build upon, unlocking additional knowledge through supplementary experimentation and data. Moreover, the model itself is flexible enough to be used as a framework for studying other biologically active metals, such as zinc or nickel, the former being another biologically important divalent ion [20]. The use of deletion strains would permit the comparison of parameter estimations to quantitatively analyze the transport and enhance the understanding of these processes.

Materials and methods
Wild-type yeast strain BY4741 were cultured in a glycerol-based growth media for twenty-eight hours. The use of glycerol ensured that the culture respired for the duration of the experiment and prevented a metabolic shift from fermentation which could alter the population and affect the outcome of the growth curve, complicating mathematical analysis. The number of cells in the culture at any given time was estimated through optical density measurements which are unitless values reported based on absorbance of a sample at a wavelength of 600 nm, known as OD 600 . Samples were diluted to ensure that they were within the linear range of the spectrophotometer for reliable measurements. The cell number factor was required to estimate the number of cells per OD 600 unit. For wild-type yeast, the average cell number factor is close to 1×10 8 cells ml −1 [21].
Isotope dilution mass spectrometry (IDMS) was used to quantify the amount of copper in a sample with high accuracy using measured isotopic compositions. A known amount of a stable 'spike' isotope is added to a known amount of sample. The mixture of spike and sample is left to equilibrate and the isotopic composition of the mixture is measured [22]. Copper has two stable isotopes, 63 Cu∼69% abundant and 65 Cu~31% abundant. Samples were weighed and then mixed with known amounts of 65 Cu-enriched spike. These mixtures were digested with concentrated nitric acid and hydrogen peroxide in a microwave to break down the biological matrix. Once in solution, the copper was isolated from the matrix by ion exchange separation [23].
The copper isotopic composition of the samples was measured on a multiple-collector inductively coupled plasma mass spectrometer (Thermoscientific Neptune MC-ICP-MS). The amount of copper was determined through the following equation, where c spike is the concentration of copper in the spike solution, m spike is the mass of the spike in the mixture, both m a,sample and m a,spike are the atomic masses of the sample and spike respectively, A 63Cu,spike and A 63Cu,sample are the abundances of 63 Cu in the spike and sample, and R spike , R sample+spike , and R sample are the isotope amount ratios.

Uncertainty assessment
Uncertainty in measurement values were evaluated through the Kragten method [24]. The Kragten method [24] is a numerical approach to assessing a standard measurement uncertainty that reduces the risk of calculation error by using spreadsheet operations to approximate the partial derivatives of the GUM method [25]. The standard measurement uncertainties are then multiplied by a coverage factor (k) of 2 to obtain an expanded measurement uncertainty (u) on a confidence interval of 95%. Major contributions to the uncertainty in the determination of the copper concentrations included the uncertainty associated with the concentration of the 65 Cu enriched solution, blank contributions from sample digestion and ion exchange, and the uncertainty of the masses of samples and spikes mixed. The copper concentration of 65 Cu enriched solution is calibrated using mixtures of copper isotopic reference materials ERM-AE633 and ERM-AE647 to obtain an uncertainty in the copper concentration of 1.2%.

Inversion algorithm
The inversion algorithm selected, fminunc, is a function contained in the MATLAB [26] programming environment for locating the minimum in unconstrained functions of multiple variables. The algorithm is based on the quasi-Newton root finding algorithm and utilizes the (BFGS) formulation [27]. The error tolerance for the algorithm was set to machine epsilon ε=2×10 -16 , the best precision for a 64-bit double-precision floating-point system.

Data and copper concentration results
Experimental copper data obtained through IDMS is presented in table 1. Variable volumes of the culture were collected in time as the size of the culture increased to ensure that samples contained numbers of cells proportional to an optical absorbance density of six. Optical absorbance data collected with a 600 nm spectrophotometer is provided as a number that may be multiplied by the cell number factor 1.0× 10 8 cells ml −1 to provide an estimate of the number of cells per volume. For an optical density of six, the number of cells in each representative sample is 6.0×10 8 cells.
Each sample was introduced to a centrifuge to pellet the cells it contained. The mass of each pellet was close to 36 mg, so this was used as a constant to correct for varying water content and multiplied by the copper concentration of each pellet (in nanograms of copper per gram of cells) to provide an absolute copper measurement (in nanograms of copper). The absolute copper content was then divided by the number of cells in the pellet, 6.0×10 8 cells, to yield the amount of copper per cell.
The growth media was assumed to be a large reservoir of copper such that the concentration was constant during the experiment. The copper concentrations from the first hour of collection were averaged to provide a mean value of 89 ng of copper per gram of growth media. During the first hour, there was 228.7 ml of growth media at a density of 1.034 g mL −1 and the concentrations were multiplied by 236.5 g to obtain the total copper. The value 21091 ng of copper was used as a constant to represent the amount of copper contained in the growth media.