Loading equine oocytes with cryoprotective agents captured with a finite element method model

Cryopreservation can be used to store equine oocytes for extended periods so that they can be used in artificial reproduction technologies at a desired time point. It requires use of cryoprotective agents (CPAs) to protect the oocytes against freezing injury. The intracellular introduction of CPAs, however, may cause irreversible osmotic damage. The response of cells exposed to CPA solutions is governed by the permeability of the cellular membrane towards water and the CPAs. In this study, a mathematical mass transport model describing the permeation of water and CPAs across an oocyte membrane was used to simulate oocyte volume responses and concomitant intracellular CPA concentrations during the exposure of oocytes to CPA solutions. The results of the analytical simulations were subsequently used to develop a phenomenological finite element method (FEM) continuum model to capture the response of oocytes exposed to CPA solutions with spatial information. FEM simulations were used to depict spatial differences in CPA concentration during CPA permeation, namely at locations near the membrane surface and towards the middle of the cell, and to capture corresponding changes in deformation and hydrostatic pressure. FEM simulations of the multiple processes occurring during CPA loading of oocytes are a valuable tool to increase our understanding of the mechanisms underlying cryopreservation outcome.

Oocytes can be used in artificial reproduction technologies. Their longevity in vitro, however, is only limited. Cryopreservation can be used to store oocytes for extended periods at ultra-low temperatures in liquid nitrogen tanks or mechanical freezers until they are needed. Both the water-to-ice phase transition as well as the drastic temperature changes during cooling and re-warming associated with cryopreservation of cells can be very damaging 1 . Membrane permeating cryoprotective agents (CPAs) like dimethyl sulfoxide (DMSO), glycerol (GLY), or ethylene glycol (EG) need to be added to mitigate the damaging effects of ice formation and the large temperature excursions 2 . Alternatively, vitrification (ice-free cryopreservation) can be used for cryostorage, which requires a combination of high cooling rates and high CPA concentrations to completely avoid ice formation 3 . For both cryopreservation approaches, the introduction of CPAs into cells needs to be done with care, because exposing cells to CPA solutions may be toxic and cause osmotic damage.
The first step in (ice-free) cryopreservation procedures for cells is to load the cells with CPAs. Cellular membranes present a barrier for CPAs to provide intracellular protection 4 . Each cell type has its own characteristic membrane permeability characteristics that define the rate at which water and solutes (i.e., CPAs) move across the membrane. The direction of water and solute transport is governed by the direction of the osmotic and solute gradient between the intra and extra-cellular environment. Transferring cells into a solution containing molar concentrations of CPAs causes cells to respond by moving water out of the cell and CPAs into the cell. Water can pass the cellular membrane much faster compared to CPAs causing water to initially move out of the cells and a concomitant reduction in cell volume. Thereafter, both CPAs and water move into the cell until equilibrium is reached between the extra-and intracellular osmolality. Cell volume reaches a new equilibrium volume that is somewhat larger as the original volume to accommodate the CPAs initially not present inside the cell. This biphasic cell volume response during the exposure of cells to CPA solutions can be captured by www.nature.com/scientificreports/ video-microscopy [5][6][7][8] . Similarly, cryopreserved cells that are initially loaded with molar concentrations of CPAs after rewarming undergo swelling during CPA removal when cells are exposed to physiological media. Different mass transport formalisms are available to fit cell volume response data to derive the characteristic cell membrane permeability to water (L p ) and the solute (P s ). Briefly, these include the so-called one-parameter (only water transport) or solute permeability model, the two-parameter model in which water and solute transport are considered as independent processes, and a three-parameter model which takes into account solute-solvent interactions during transport across the cellular membrane 9,10 . Once cell-specific L p and P s values are known, they can be used to predict volume responses during CPA loading and unloading at different CPA concentrations and, if the activation energies of water and CPA permeability are known, at different temperatures 11 .
The analytical mass transport models mentioned above do not take CPA diffusion inside the cell into account and do not provide spatial information of the complex processes associated with CPA (un)loading of oocytes. A few recently developed mathematical models, however, have implemented spatial information of intracellular CPA concentration 12,13 . The finite element method (FEM) can be used to develop phenomenological simulations of mass transport and mechanical behavior in a 3D framework 14,15 . The FEM has been used to model the response of oocytes exposed to mechanical stress in a micromanipulator set up 16 , and to simulate mechanical properties of human oocytes 17 . The FEM subdivides a system into smaller parts (the finite elements) using spatial discretization leading to the emergence of a so-called mesh that represents the finite element problem. The numerical treatment is developed based on a single element and then, it is extended to the entire discretized domain through the standard procedure of assembly based on the element connectivity.
The main aim of this study was to develop a phenomenological FEM model that can be used to design CPA loading protocols for equine oocytes. Parameters of the oocyte membrane permeability to water and EG were inferred from previously obtained video microscopic imaging results 18 . First, the two-parameter formalism was used to model the volume response and concomitant intracellular CPA concentration during the exposure of oocytes to EG solutions of varying concentrations and at different temperatures. The FEM was used to model deformation, mechanical stress, and intracellular CPA concentrations at different locations in the cell.

Results
Volume response and CPA uptake of oocytes during exposure to CPA solutions. Figure 1A shows the volume response of equine oocytes exposed to 1.5 mol L −1 EG at 22 °C. Cell volume response data were fitted using Eqs. (1) and (2) and the parameters listed in Table 1 to derive L p and P s values as previously described 18 . The membrane permeability values for water and EG and their concomitant activation energies (see Table 1) were used to simulate oocyte responses during exposure to varying EG concentrations as well as varying temperatures ( Fig. 1B−E).
Oocytes exhibit a biphasic cell volume behavior when exposed to solutions containing molar concentrations of EG. First, the volume decreases, due to water moving out of the cell, where after oocytes return to their original volume due to water and solutes/EG moving into the cell. Both the initial rate of cellular dehydration as well as the minimum volume that is attained during EG exposure appears to increase with increasing EG concentrations (Fig. 1B). The minimum volume that is attained after approximately 30 s of EG exposure should not exceed the osmotic tolerance limits of the cell. The intracellular EG concentration initially shows a rapid increase coinciding with the volume decrease during EG exposure where after the intracellular EG concentration gradually approaches the extracellular EG concentration (Fig. 1C).
Temperature has a profound effect on the cell volume response and EG uptake. Both the initial rate of cell dehydration during EG exposure and the rate at which the cell returns to its original volume decrease with decreasing temperature (Fig. 1D). Moreover, the maximum extent of cell dehydration during EG exposure (i.e., the minimum volume that is attained during EG exposure) is greater at lower temperatures. The slower volume response at lower temperature coincides with a slower EG uptake and increases the time it takes to reach osmotic equilibrium (Fig. 1E).
Use of FEM simulations to model the response of oocytes during exposure to CPA solutions. Oocyte deformation (i.e., volume changes) and associated changes in hydrostatic pressure during EG exposure were modeled using FEM simulations. Figure 2A-G illustrates different views of an oocyte held with a micro-capillary as obtained with FEM analysis, and the flow chart that was used for the simulations (details are described in the 'Materials and methods' section). FEM simulations were done using the parameters listed in Table 2, which were optimized using the two-parameter formalism. Simulations were done for different extracellular EG concentrations to visualize cellular deformation and hydrostatic pressure in the mesh elements reflecting different cellular locations ( Fig. 3A−C). The FEM simulations show a biphasic cell deformation (i.e., volume) response at different EG concentrations, which closely resemble the volume responses obtained with the twoparameter model (indicated as dotted lines in Fig. 3D). The hydrostatic pressure also shows a biphasic response; the pressure first increases to reach a concentration-dependent maximum after which the pressure gradually decreases again (Fig. 3E). The hydrostatic pressure profiles precede the deformation profiles; the maximum pressure is attained after about 15 s, whereas the minimum cell volume is attained after approximately 30 s. Figure 4A shows FEM simulations of the spatial EG distribution in an oocyte during EG exposure. At the end of the simulation period at 4 min, the intracellular EG concentration near the membrane surface is clearly greater compared to the layers towards the center of the oocyte (Fig. 4B). Figure 4C shows that the average intracellular solute concentration during EG exposure (obtained by integrating the EG concentrations in the mesh elements depicting the oocyte) increases substantially slower versus time compared to the profile obtained using the two-parameter model (indicated as dotted lines). The two-parameter formalism only takes diffusion over the plasma membrane into account and not the EG diffusion in the cell. This explains why the kinetics of www.nature.com/scientificreports/ EG uptake near the membrane surface resembles the kinetics obtained by the two-parameter formalism, whereas the EG concentration versus time profiles in the layers towards the center of the cell show a much slower and more gradual increase in EG concentration over time (Fig. 4D). After thawing, cells are exposed to physiological solution (in vitrification procedures this is done in a stepwise fashion). The CPA that is still inside the cells result in a concentration gradient and an osmotic imbalance resulting in CPA unloading. Just like the CPA uptake, CPA unloading is also characterized by a biphasic cell volume behavior. The presence of intracellular CPAs after thawing causes an influx of water and diffusion of CPAs to the extracellular buffer. First, the volume increases, due to water moving into the cell, where after oocytes return to their original volume due to water and solutes/EG moving out of the cell (Fig. 5A). The intracellular EG concentration initially shows a rapid decrease coinciding with the volume increase during isotonic solution exposure where after the intracellular EG concentration gradually approaches zero (Fig. 5B). Just like as has been observed with CPA uptake (Fig. 4C), the CPA efflux curves obtained with the FEM simulations are delayed compared to those obtained with the two-parameter model.
It is known that symmetric systems can be treated using 2D or even 1D approximations depending on the extent of symmetricity. However, this simplification is not applicable to non-symmetric cases and one needs to do a full 3D analysis in order to obtain the distribution of the variables in the entire domain. For instance, if a cell is subjected to non-symmetric boundary conditions, one expects to get a non-symmetric response. This can be simulated using non-symmetric boundary conditions of the CPA concentration on the exterior surface of the cell (Fig. 6A). The boundary conditions for this test case are essentially similar to those of the symmetric case shown in Fig. 2A except for the solute concentration at the exterior surface. One can see that part of the cell (shown in light green in Fig. 6A) is exposed to a lower concentration (0.5 C ref ) while the rest of the cell surface is exposed to the maximum reference concentration C ref (shown in dark green in Fig. 6A). From a physical point of view, this example mimics a situation in which the permeability of the cell membrane is weaker or stronger in some areas leading to a non-uniform flux of both CPA and water into and out of the cell. Particularly with Table 1. Parameters and values used for 2-parameter fitting and modeling of equine oocyte volume responses towards solutions containing various EG concentrations. w: water, s: permeating solute, n: non-permeating solute, o: initial value. e: extracellular, i: intracellular. *L p and P s values were experimentally derived for exposure to 1.5 M EG, at room temperature, as described previously 18 . These values were also applied for calculating volume changes in response to exposure to other EG concentrations. **The values were obtained from the literature 8 . ***Furthermore, these values were used in Eqs. (7) and (8)  www.nature.com/scientificreports/ glycerol, which has a relatively low P s , non-uniform shrinkage has been observed (Fig. 6B). The FEM simulation (Fig. 6C) shows a qualitative comparison with the microscopic observations.

Discussion
Analytical mass transport models have been widely used to simulate cell volume responses and CPA uptake during the exposure of cells to CPA solutions 5,6,8,19 . Modeling predictions can be used to avoid conditions that cause cells to shrink or swell beyond their osmotic tolerance limits and to compute the time needed to reach a sufficiently high enough intracellular CPA concentration to protect against freezing injury 11 . For vitrification approaches modeling can be used to estimate the time needed to reach maximal dehydration of the cell prior to www.nature.com/scientificreports/ cooling, while keeping the exposure time as short as possible to reduce toxicity effects 20 . Precise timing is of the essence to maximize the cryopreservation outcome. Exposing cells to anisotonic solutions above or below their tolerance limits results in a rapid decline in cell viability. The osmotic tolerance limits define the maximum tolerable volume excursions that a cell can withstand via the Boyle van't Hoff equation. For red blood cells, cell volume excursions should be maintained below 1.7 times the isotonic osmotically active volume to avoid hemolysis 21 .
Oocyte deformation is a major cause of potential cell damage which can lead to disruption of the spindle apparatus and failure in embryo development after in vitro fertilization 22,23 . Therefore, for oocytes, the osmotic tolerance limits are not only defined by the apparent membrane integrity after exposure to anisotonic solutions, but also by their developmental potential. Porcine and equine oocytes respond as linear osmometers when they are exposed to anisotonic solutions ranging from, respectively, 180 to 800 mOsm kg −1 and 200 to 600 mOsm kg −1 suggesting that they can tolerate exposure to anisotonic solutions in these osmotic ranges, i.e., the oocyte membrane remains intact 4,18 . However, damage in the oocyte spindle apparatus only becomes apparent in its developmental potential, i.e., the ability of an oocyte to form blastocysts after in vitro fertilization or intracytoplasmic sperm injection. A study on porcine oocytes showed that only about half of the oocytes maintain their developmental potential after exposure to anisotonic solutions only slightly above or below isotonic conditions 24 .
In the 2-parameter formalism, only the plasma membrane (i.e., outer membrane) is included as diffusion barrier. As a consequence this model yields a rapid increase of the intracellular CPA concentration reaching a plateau within 30−60 s. The intracellular environment is not solely comprised of fluid, but also includes organelles all surrounded by their own membranes (i.e., apolar diffusion barriers). Use of the FEM model actually may give a more realistic estimation of the intracellular CPA concentration during CPA loading, namely a slower increase in intracellular CPA concentration. For cryopreservation purposes, the mechanical response of the cell is particularly important, i.e., the time it takes for the cell volume to attain equilibrium again. The simulations show that the mechanical response of oocytes during CPA exposure (shrinkage followed by swelling) actually is not yet fully completed in 4 min.
The two-parameter formalism only requires knowledge of the membrane permeability characteristics and dimensions of the cell. However, exposing a cell to osmotic imbalance causes a complex series of events including deformation as well as mass transport in the extracellular space, across the cell membrane and within the intracellular spaces. The FEM was used here to develop a phenomenological continuum model to simulate strain-stress behavior of an oocyte kept in position by a holding pipette during exposure to a CPA solution without using analytical solutions and prior knowledge of molecular and microscopic structures. The parameters used in the FEM model were optimized using the two-parameter formalism so that the characteristic biphasic cell volume behavior of both models resembled each other, even though the volume responses were not identical. Simulations of the EG concentrations in the mesh elements depicting the oocyte provided insights in the kinetics of the CPA loading process in different cellular locations of the cell. Close to the membrane surface, the uptake kinetics resembled the kinetics obtained using the two-parameter formalism, whereas the uptake in locations closer to the center of the cell, show a delayed and slower CPA uptake kinetics different as predicted by the two-parameter formalism. The FEM model was developed for EG, but can easily be modified for other CPAs by changing the parameters in the FEM algorithm, particularly the diffusivity and the initial expansion or contraction coefficient.
Summarizing, whereas analytical transport models can be used to simulate average oocyte volume responses and concomitant intracellular CPA concentrations during exposure of oocytes to CPA solutions, FEM models provide more detailed and complementary information of the CPA loading process. FEM simulations can be used to depict spatial differences in CPA concentration during CPA permeation and to capture corresponding changes in deformation and hydrostatic pressure. Modeling the response of oocytes exposed to CPA solutions can be used to define suitable conditions to safely load oocytes with CPAs, and to remove the CPAs again after thawing.

Conclusion
The objective of this 3D continuum model is to perform a full analysis of the oocyte response in the presence of a CPA solution. One can extend the model and incorporate thermal effects on intracellular CPA distribution during freezing or thawing. Spatial information on local CPA concentrations can in turn be used to predict the likelihood of intracellular vitrification and/or intracellular ice-formation. Such models are particularly useful if they are supported and verified by microscopy results. Once the model is validated, it can be used to conduct www.nature.com/scientificreports/ www.nature.com/scientificreports/ in silico experiments that can be used to optimize cryopreservation or vitrification protocols for cells, which is substantially cheaper than in vitro investigations. In our upcoming work, we plan to incorporate local structures including the plasma membrane and organelles in conjunction with the diffusion process within the bulk of the cell. The explicit modelling of water content using 'porous media framework' as well as 'thermal effects' can also be introduced into the model.

Materials and methods
Experimental assessment of oocyte membrane permeability parameters. Data from a previous study 18 were used to develop analytical and FEM models that can be used to simulate oocyte volume responses during exposure to CPA solutions. In this previous study 18 , equine oocyte volume responses were analyzed upon perfusion with saline supplemented with 1.5 mol L −1 permeating cryoprotective agent (i.e., ethylene glycol/EG www.nature.com/scientificreports/ or glycerol/GLY). Volume responses were recorded and analyzed, using an inverted microscope with a micromanipulator setup. Briefly, oocytes were transferred to a 10 μL droplet of isotonic medium on a culture dish, which was placed on the microscope-micromanipulator stage. A holding capillary was used with enough negative pressure to hold the oocyte at a fixed position. A syringe was used for adding 1.5 mL CPA solution onto the droplet in which the oocyte was held. The two-parameter transport formalism 10 was used to fit oocyte volume versus time plots to derive the membrane permeability to water (L p ) and for EG (P s ) as described in our previous paper 18 . Fitting was done using MATLAB software (Mathworks, Natick, MA, USA). The cell volume (Vc) change during the perfusion process, at a particular temperature (T), can be described as a function of the water and solute volume (V w and V s , respectively): furthermore: where, V s is the partial molar volume of the CPA, M the osmolality, w refers to water, s to permeating and n to non-permeating solute, 0 refers to the initial value, while e and i refer to external and internal cellular locations, respectively. R represents the universal gas constant and cell specific parameters include the area (A), isotonic volume (V 0 ), the osmotically inactive volume (V b ), and osmotically active volume ((1 − V bf ) × V 0 ) with V bf referring to the fractional value of cell solids. It was assumed that CPA solutions behave as ideal solutions.
The temperature dependence of the transport parameters can be described using the following Arrhenius relationship: www.nature.com/scientificreports/ www.nature.com/scientificreports/ where, rates and temperatures are given per s and in K, respectively, while L pg and P sg are the reference values of L pg and P sg at the reference temperature T 0 (i.e., 273.15 K which equals 0 °C), and E Lp and E Ps represent the activation energies for water and solute (CPA) transport across cell membrane, respectively.

Modeling oocyte volume responses and intracellular CPA concentrations during exposure to different CPA concentrations and temperatures.
Membrane transport parameters, L p and P s , of equine oocytes were previously studied using the oocyte volume response upon exposure to 1.5 mol L −1 EG 18 . The obtained L p and P s parameters were used to simulate oocyte volume responses and concomitant intracellular CPA concentrations upon exposure to varying EG concentrations and at different temperatures. A detailed listing of the parameters and values that were used for modeling equine oocyte volume responses and concomitant intracellular solute concentrations according to the two equations of the two-parameter formalism is presented in Table 1. A user-friendly interface was created with the App-designer included in MATLAB. The code includes a section for calculating the L p and P s values at varying temperature (according to Eqs. 7 and 8): furthermore, it is linked to an Excel data sheet from which the following parameters are imported: L pg = cell membrane hydraulic permeability (L p ) value at 273.15 K. P sg = cell membrane solute permeability (P s ) value at 273.15 K. E L P = activation energy of cell membrane hydraulic permeability. E P s = activation energy of cell membrane solute permeability. R = universal gas constant. T = room temperature. T 0 = initial temperature (0 °C = 273.15 K). In a separate section, the oocyte volume responses versus time are then calculated using: where: V w0 = isotonic (initial) water volume. V b = osmotically inactive cell volume. where: M i s = internal permeating solute osmolality, V n = normalized volume. Finally, the numeric data of the calculated cell volume and solute concentration versus time plots are exported to an Excel data sheet, and plots were created using Sigmaplot for Windows version 13.0 (Systat Sotware Inc., http:// www. systat. de/ index. html).
Modeling of diffusion-induced deformation of oocytes during exposure to different CPA concentrations using ANSYS. Using ANSYS (ANSYS Workbench, academic license version 19.2, https:// www. ansys. com/), as a solver based on FEM, simulations were conducted in which cell volume change is considered to be the result of diffusion-induced mechanical deformation. The cell volume change, from a physical point of view, originates from variations in the cell water content which is dictated by the cell membrane permeation and osmotic pressure. In our model, we only took into account the diffusion of the solute, whereas the osmoti- www.nature.com/scientificreports/ cally driven water flow was not explicitly incorporated in the model. Nevertheless, the indirect and implicit impact of the volume change due to the cellular water content was introduced into the deformation process. It is a phenomenological model in which one can interpret the volume change of the cell as a process depending on the concentration of the solute, namely CPA. As a matter of fact, the dependency of water content, as the one that is directly responsible for the volume change, on the CPA concentration plays the key role. If we were to account for the water 'explicitly' , we would have to extend the model to a poro-elastic one in the framework of porous media theory. Hence, for the sake of simplicity, we exclude this part by assuming a phenomenological model and directly connecting the CPA concentration to the local volume change (mechanical strain). The coupling between the deformation and the diffusion-driven expansion was feasible through developing a computer code using ANSYS Parametric Design Language (APDL) 25 . The applied model is transient, meaning that variables are dependent on not only the position but also the time. The inertial effects (i.e., mechanical acceleration), however, were not taken into account due to the fact that the process is relatively slow compared to systems without an internal damping effect (e.g., a load in static structural analyses does not change or changes slightly that do not disrupt the steady-state conditions) and hence dynamics does not play a role here. Unlike mechanical processes, the solute diffusion phenomenon is time-dependent.
Appropriate boundary conditions were chosen to mimic the experimental set up. Two field variables are considered, namely mechanical deformation and CPA concentration. The oocyte is held at a specific region using a capillary with a defined diameter, while the solution containing specific solutes/CPAs is added ( Fig. 2A). From a mathematical point of view, the oocyte is mechanically subjected to so-called Dirichlet boundary conditions. It means that the deformation components are set to zero at the zone kept by the capillary (indicated in red in Fig. 2A). The rest of the cell is free to undergo deformation due to swelling/shrinkage. Similar to the mechanical field, the boundary condition of the diffusion equation is of Dirichlet type at the regions submerged in the CPA solution. It means that the maximum solute concentration (denoted by C ref ) is prescribed on the exterior surface of the cell (shown in dark green in Fig. 2A). Nevertheless, in the cellular region that is in contact with the capillary (shown in red in Fig. 2A), a zero flux condition is enforced. In mathematical terms, it corresponds to a Neumann type boundary condition in which the concentration gradient is prescribed.
The simulation procedure can be divided into three sections, namely: pre-processing, solving and postprocessing modules. In the pre-processing part, the following inputs are defined: the initial conditions, types of elements used, diffusion coefficients, saturation and reference concentrations, the elasticity module and Poisson's ratio, and diffusion expansion coefficients. The solid227 element was chosen to perform combined structural diffusion analysis. When used in coupled-field analyses having structural and diffusion DOFs; this element has ten nodes with up to six degrees of freedom per node and supports the effects of diffusion strain and hydrostatic stress migration (i.e., transport of particles due to a hydrostatic stress gradient).
In the solving-section of the program, first the boundary conditions are defined and it will be determined how long the simulation will take. For being able to compare the simulation with the experimental data (i.e., oocytes were kept in place with a holding capillary), some points of the sphere were fixed resulting in a restricted shape change. A summary of the mathematical framework behind the model is described hereafter. As a convention, unlike scalar variables, vectorial variables are represented in bold, while the matrices are distinguished using squared brackets.
The mechanical strain is defined using the gradient of displacement vector (denoted by u) as follows: in which the Nabla operator (∇) signifies the spatial gradient. Assuming a simple Fickian transport, the mass flux vector J for the diffusion of the solute is described by: www.nature.com/scientificreports/ where J represents the diffusion flux density and [D] as the diffusivity matrix is defined according to: in which D XX , D YY , and D zz are diffusivity coefficients in the X, Y, and Z directions, respectively. Similar to the diffusion expansion coefficient, the diffusivity coefficients are taken to be identical and consequently it leads to an isotropic diffusive behavior, namely D XX =D YY = D ZZ = D.
It is assumed that the diffusion expansion parameter β is a function of the solute concentration which controls the water content. It is obvious that the local volume of a material point changes due to a change in the volume fraction of the water there. Since the osmotic flow of the water is not explicitly present in our model, the diffusion expansion parameter is directly assumed to be a function of the solute concentration. For the sake of simplicity, this dependency is assumed to be of the exponential type according to: where β 0 and C 0 describe the dependence of β versus dimensionless concentration C.
It should be noted that C is a dimensionless representation of the solute concentration according to: Figure 2B illustrates the variation of β as a function of C.
One should notice that Eq. (21) is only a choice that contains two meaningful parameters: Initial value ( β 0 ) and decay rate (C 0 ). Here, the water content variation as the main driver of swelling/shrinkage is translated into the evolution of parameter β. From a physical point of view, the water content depends on the availability of CPA concentration, namely β = β(C). Since the CPA concentration at any point is a function of time due to the transient diffusion equation, namely C = C(x,t), it implies that the parameter β is a function of time. The overall volume of the cell changes in the course of time. The volume simulations obtained using the 2-parameter formalism form the basis to develop the FEM model. Hence, β 0 and C 0 have been calibrated using the volume curves obtained with the 2-parameter formalism, which in turn are based on experimentally determined volume response curves of oocytes exposed to a 1.5 mol L −1 EG solution 18 (i.e., relative cell volume change versus time plots in Fig. 3D for 1.5 mol L −1 ). In other words, the optimal value for these two parameters are found in such a way that fit experimental observations of the volume change over time. Once β 0 and C 0 have been estimated, one can use them similar to other invariant parameters (E,ν,D) in order to predict the behavior of the system for other conditions.
The balance of mass for the solute concentration C is described by the following diffusion equation: where the overhead dot stands for time derivative and the symbol ∇. is defined as the divergence operator. In addition to the balance of mass, one should take into account the equilibrium of mechanical forces that is governed by: where σ signifies the mechanical stress tensor defined in Eq. (17). It is apparent that, unlike Eq. (23), the mechanical governing Eq. (24) is assumed to be time-independent. Applying the variational principle to the structural (Eq. 24) and diffusion equations (Eq. 23) coupled via the constitutive equations (Eqs. 17 and 19), one can obtain the following finite element matrix equations for the structural-diffusion analysis:   www.nature.com/scientificreports/ The nodal displacement vector u and nodal concentration C are the unknowns in the finite element problem in hand. On the other hand, the sum of the element nodal force F and nodal diffusion flow rate vector R, which constitute the right hand side of the finite element formulation, are known due to the prescribed boundary conditions. Figure 2C-F illustrate the different views that can be obtained with finite element analysis in ANSYS, namely: cross sections at different locations (Fig. 2C,D), meshed shapes (Fig. 2E), and the holding point area (Fig. 2F). In the post-processing section of the program, one can generate graphical outputs of the simulation such as contours and animations. The ANSYS workflow is illustrated in Fig. 2G. The ANSYS model is based on 5 parameters: E, ν, D, β 0 , C 0 reported in Table 2. Because experimentally captured oocyte volume excursions took place during a 4 min time interval 18 , we kept this time for the ANSYS simulations for concurrent evaluation of cell deformation and CPA concentration.