Predicting the survival rate of mouse embryonic stem ... - SAGE Journals

2 downloads 0 Views 1MB Size Report
Predicting the survival rate of mouse embryonic stem cells cryopreserved in alginate beads. S Sambu, X Xu, H Ye, and Z F Cui*. Institute of Biomedical ...
1092

Predicting the survival rate of mouse embryonic stem cells cryopreserved in alginate beads S Sambu, X Xu, H Ye, and Z F Cui* Institute of Biomedical Engineering, Department of Engineering Science, University of Oxford, UK The manuscript was received on 12 October 2010 and was accepted after revision for publication on 8 June 2011. DOI: 10.1177/0954411911418568

Abstract: Stem cell cryopreservation in three-dimensional (3D) scaffolds may offer better protection to cells leading to higher survival rates. However, it introduces heterogeneity in cryoprotective agent (CPA) concentrations, durations of exposure to CPA, and freezing and thawing rate within constructs. This paper applies a mathematical model which couples the mass transport of dimethyl sulphoxide (DMSO) in a cell-seeded spherical construct and cell membrane transport into mouse embryonic stem cells (mESCs) to predict overall cell survival rate (CSR) based on CPA equilibrium exposure times (tE) and concentrations. The effect of freeze-concentration is also considered. To enable such a prediction, a contour plot was constructed using experimental data obtained in cryopreservation of cell suspensions with DMSO at a cooling rate of 1 °C/min. Thereafter, the diffusion in the alginate bead and the membrane transport of CPA was numerically simulated. Results were mapped onto the survival rate contours yielding ‘predicted’ CSR. The effects of loading time, hindrance, construct radius, and CPA concentration on predicted CSR were examined. From these results, an operation window with upper and lower tE of 12–19 min (for 0.6 mm radius beads and 1.4 M DMSO) yielded an overall viability of 60 per cent. The model predictions and the best experimental cryopreservation results with encapsulated mESCs were in agreement. Hence, optimization based on postthaw CSR can accelerate the identification of cryopreservation protocols and parameters for maximizing cell survival. Keywords: cryopreservation, mouse embryonic stem cells, alginate encapsulation, mass transfer, cell membrane transport, DMSO, operation window, freeze-concentration

1

INTRODUCTION

Stem cell cryopreservation with quantitatively predictable outcomes may ensure an adequate and timely supply of cells for regenerative medicine. Cryopreservation combines the cryogenic arrest of biochemical reactions and the freeze-induced lack of water to hold cell biology in abeyance. There are two main cryopreservation methods: slow freezing and vitrification. Slow freezing uses a controlled cooling rate and low cryoprotective agent (CPA) *Corresponding author: Institute of Biomedical Engineering, Department of Engineering Science, University of Oxford, Oxford, UK. email: [email protected]

Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

concentrations; vitrification uses a rapid cooling rate and a high CPA concentration. Although slow cooling and vitrification of cells is often done in suspension, developments in tissue engineering have motivated the cryopreservation of cell–material constructs [1, 2]. Additionally, because the preservation of viability and particular functions of stem cells is a prerequisite for their clinical and industrial applications, new techniques are needed to achieve those aims simultaneously. Biomaterials may be a way to achieve these two aims. For example, they can minimize the effects of ice crystal growth during cryopreservation [3]. Additionally, by providing a substratum, three-dimensional (3D) cryopreservation can regulate the changes in cell volume during cryopreservation. The stabilization of

Predicting the survival rate of mouse embryonic stem cells

cell phenotype by using biomaterials has been demonstrated during culture [4]; such biophysical cues may be extended to cells at the cryopreservation stages before and after freezing, leading to low levels of post-thaw spontaneous differentiation. One useful polymer for cryopreservation is alginate. Alginate is a block co-polymer of guluronic and manuronic acid. It is an attractive biomaterial because of its low immunogenicity, easy sol–gel transition, abundance, low cost, and an amenable chemistry allowing for the design of specific cell– material interactions. Its rapid sol–gel transition enables the extraction of cells using chelators for the divalent ions. To provide the 3D microenvironment for cryopreservation, mouse embryonic stem cells (mESCs) can be encapsulated in alginate gel beads. In encapsulation, each cell is surrounded by a polymer gel. Specifically, spherical drops of cells suspended in alginate can be introduced into a calcium chloride solution; calcium ions will cross-link the alginate allowing the formation of calcium– alginate beads with entrapped cells. One of the earliest uses of encapsulation was the culture and immunoisolation of pancreatic islet cells for the production of insulin [5]. Since then, various cell types have been used and the necessity for longterm storage has motivated investigations of 3D cryopreservation techniques. To minimize the ice crystal size during cooling, a slow cooling rate was achieved using a programmable freezer [6]. Although vitrification of encapsulated cells has been attempted [7], it has proven cumbersome due to the high viscosity of the cryopreservative [8] and the difficulty in achieving uniform ultra-fast freezing. These limitations of vitrification are eliminated when slow cooling is adopted, leaving the optimization of the loading conditions as a final barrier to achieving improved cryopreservation results [9]. Models have been developed to better understand the underlying physical phenomena. Most models are based on the mass transport of CPA and water for engineered tissues [10] and for specific tissue explants [11]. While these models describe the detailed physical aspects of cryopreservation, the impact of physical phenomena on postthaw cell survival rate (CSR) was not directly demonstrated. In this paper, an attempt has been made to predict the CSR after the cryopreservation of cells encapsulated in alginate. Dimethyl sulphoxide (DMSO) was the chosen CPA and mESCs were the model cells. The survival rate was modelled based on DMSO concentration and equilibrium exposure time; comparisons were made to models explicitly

1093

incorporating freeze-concentration showing a similarity in survival rate predictions. The criteria were derived using experimental data from the cryopreservation of cell suspension using slow-cooling methods. Both parameters can be predicted using a mass transfer model and hence the overall CSR for the whole bead can be evaluated. 2

METHODS FOR THE IMPLEMENTATION OF THE NUMERICAL MODELS

2.1 Procedures for the cryopreservation of encapsulated cells Cell encapsulation was performed in three major steps: (a) cells were harvested in suspension; (b) cells were then resuspended in alginate at a given cell density, r; (c) finally, the alginate–cell mixture was ejected (drop-wise) into a calcium chloride solution for a preset period of time to form spherical beads and the beads were washed. The encapsulated cells were then cryopreserved and thawed via a five-step process: (i) the beads were immersed in a CPA solution at a predetermined CPA concentration, duration of exposure, and temperature; (ii) the temperature was decreased at 1 °C/min to 280 °C (at this stage, the beads are in contact with the bottom and sides of the vial); (iii) the beads were then stored at 280 °C; (iv) the beads were thawed in a 37 °C water bath for 1 min; (v) CPA was gradually replaced with pre-warmed culture medium over a 10 min period. 2.2 Implementation of the diffusion model The objective of this model is to predict the CPA distribution within the beads, the intracellular CPA transient, and eventually link this distribution to the post-thaw CSR. Parameters and variables used in the development of this model are provided in Table 1. The focus of this study is on ‘CPA loading’ (step i) and its impact on post-thaw CSR. CPA loading begins at t = 0, when the bead first makes contact with the CPA solution. CPA diffuses into the bead, reaches the cells, crosses the cell membrane and protects the cells during the cooling stage (step ii). CPA loading stops when cooling begins; the duration of CPA loading is termed the ‘loading time’. Cell viability depends on the local CPA concentration (C(r, t)) and the ‘equilibrium exposure time’ (tE). For 3D cryopreservation, tE refers to the Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

1094

S Sambu, X Xu, H Ye, and Z F Cui

Table 1 Reference values for selected parameters Parameter

Full name

Reference values

References

D (DMSO) D/De (alginate) j tL Ti Lp ELP Pd ref Ea D (DMSO) j tL T (t = 0) Lp Vw0 Vc nbvw Is; Ishet O k

Diffusivity Hindrance Construct radius Loading time Loading temperature Hydraulic conductivity Activation energy for hydraulic conductivity Solute permeability Activation energy for solute permeability Diffusivity Construct radius Loading time Loading temperature Hydraulic conductivity at 295 K Initial volume of cell water Instantaneous volume of the cell (c) Volume fraction of water in extracellular space Ice nucleation rates Kinetic constant for the nucleation rate Thermodynamic constant for the homogeneous nucleation rate Thermodynamic constant for the heterogeneous nucleation rate Thermal conductivity of the seeded capsule Latent heat of fusion Melting coefficient Initial concentration of permeating solute Initial non-permeating solute concentration CPA loading time Isosmotic cell volume Solid fraction in mouse ES cell Cell density Volume-normalized heat capacity

4.89 3 10210 m2/s at 4 °C 1.2 0.6 mm 1800 s 277.15 K 6.83 3 10214 m/Pa per s at 4 °C 59789 J/mol 7.65 3 1028 m/s at 4 °C 64726 J/mol 4.89E 3 10210 m2/s at 4 °C 0.6 3 1023 m 1800 s 277.15 K 6.83 3 10214 m/Pa/s 349.6 um3 (50.3 per cent of cell volume) 695 um3 Variable Variable 23.7 3 1029 m2/s 3 3 10210 K5

[12, 13] [13, 14]

khet K L c Msi Mni tL Viso Vb r rm Cp

time that cells, at a given radius of the alginate bead, are exposed to CPA at a concentration equivalent to the CPA concentration at the surface of the alginate bead (‘surface concentration’, Cs). For cells in suspension, tE is defined as the time that cells are exposed to the cryopreservative; for cells cryopreserved in suspension, tE and the loading time are equal. Having identified the key parameters and variables involved in the model, the development of the model can proceed. The cell-seeded alginate construct is immersed in a cryopreservation solution at Ti; the CPA diffuses into the construct with a diffusivity De. At the end of CPA loading, the samples are slow cooled. The basic assumptions of this model are as follows 1. The alginate construct is modelled as an isotropic, inert, and biocompatible sphere. 2. The only mechanism influencing the CPA concentration distribution within the construct during CPA loading is diffusion and membrane transport into cells. 3. The impact of freeze-concentration during slow cooling is implicitly accounted for by the splining model that relates loading history to Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

4.6 3 1029 K5 1.156

2.24 + 5.975 3 10 *(2732T) 330 3 106 J/m3 1 0 Osmoles/kg 0.29 Osmoles/kg 1800 s 695 3 10218 m3 0.497 1 3 106 cells/ml 43 106 (liquid phase); (7.16T + 138) 3 103 (solid); 1.02 3 106 (organic matter) J/m3 K

5.

6.

7.

[17] [17] [18]

23

4.

[15] [15] [15] [15] [12, 13] – – – [15] [15] [15, 16] [15, 16]

[10] [10] [10] – [15] – [15] [15] – [10]

post-thaw cell survival rate. (A comparison will be made between this model and an alternative model that explicitly accounts for freezeconcentration.) At any time, the concentration of CPA does not change with either of angles u or u of the spherical coordinate system; diffusion only occurs radially. Individual cells are evenly dispersed within the alginate construct and intracellular concentration is modelled using a two-parameter model as described by Kleinhans [19]. The cryopreservation solution is well stirred prior to CPA loading. The volume of the constructs is insignificant compared to the volume of the cryopreservation solution. Hence, the solution is homogeneous and the addition of cell-seeded constructs to the solution does not alter the concentration of the CPA in solution, respectively. The resistance to mass transfer in the cryopreservation solution is negligible. Therefore, the DMSO concentration at the surface of the construct is constant and equal to the concentration in the surrounding solution.

Predicting the survival rate of mouse embryonic stem cells

1095

8. The diffusion coefficient of the CPA is only affected by loading temperature and hindrance of the biomaterial but not by the dynamic viscosity of the solvent.

all simulations, the time step was kept above 2.3 3 10212. Discretization errors, if any, are reported during run time. The numerical solution was validated by comparison to an analytical solution.

2.3 Mathematical development of the diffusion model

2.5 Implementation and validation of the membrane transport model

The CPA loading can be described by the mathematical equations (1) to (4). The governing equation is Fick’s second law in spherical coordinates (note D/H is equal to De the effective diffusivity of CPA in the construct; the hindrance, H, increases with increase in cell density)

Cell membrane transport is then modelled onto the diffusion model to obtain the intracellular CPA concentration. Equations (5) to (9) were previously described by Kleinhans [19]. The cell volume has three chambers – water volume (Vw), solute (DMSO) volume (Vd) and ‘cell solid and non-osmotic’ volume (Vb).

  ∂c 1 D ∂ 2 ∂c = r ∂t H r 2 ∂r ∂r

(1) Vc = Vw + Vd + Vb

The initial conditions are C ð0, r Þjr\R = 0; C ð0, jÞ = Cs

The solute volume can be solved as follows (2)

The following boundary condition (at r = j) was applied C ðt, jÞ = Cs

(3)

There is no concentration gradient at the centre of the construct because of the symmetry of the problem. ∂C j =0 ∂r r = 0

(5)

(4)

2.4 Numerical solutions and validation The numerical solution was implemented in a Matlab (R2008a)TM environment. A PDE solver, pdepe, is used to discretize, fulfil initial-boundary conditions, and integrate the discretized equations with respect to time. Pdepe is a function provided in Matlab (i.e. not part of a toolbox) that solves initialboundary value problems with at least one parabolic equation. The discretization formula for the method of lines solution has previously been described [20] – briefly, it uses a Galerkin/Petrov–Galerkin method to spatially discretize the problem. Ode15s is used to perform the time-integration because it can change the time step and numerical integration formula dynamically to suit the nature of the problem. It is hence suitable should stiff problems arise [21, 22]. All the simulations in this report were performed with relative and absolute tolerance at 1 3 1023 and 1 3 1026 respectively. Errors were hence kept to the minimum of either 0.1 per cent or 1 3 1026. In

∂Vd  ∂Nd =Vd ∂t ∂t

(6)

The solute flux equation will then become coupled onto the diffusion equation as follows (initial and boundary conditions are unchanged)     ∂C 1 D ∂ 1 ∂Nd 2 ∂C = r  ∂t H r 2 ∂r ∂r Vw ∂t The solute in this model is the CPA. The two models can be linked as follows C(r, t) = Mde

(7)

The osmotic pressure difference across the semipermeable cell membrane drives the water flux.    ∂VW Nd VW0 i =  Lp ART M e   Mn0 ∂t VW VW

(8)

The solute flux is again driven by the osmotic pressure difference between solutes on either side of the semi-permeable cell membrane.   ∂Nd Nd = Pd A Mde  ∂t VW

(9)

The coupled equations for the membrane transport model were integrated using the Ode15s solver. The membrane transport model was validated by comparison to experimental data for mouse C57BL/ 6 ES cell line [15]. Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

1096

S Sambu, X Xu, H Ye, and Z F Cui

2.6 Mathematical development of the heat transfer model to account for freeze-concentration and post-thaw CSR The heat transfer model is based on the conduction equation and includes a term to account for ice formation. The model also takes into account the physical and cell-mediated mass transfers that proceed in the liquid phase prior to full solidification via mass balances. The aim of the model is to link temperature transients to freeze-concentration (FC); exposure is measured by calculating the area under the molarity–time graph and (with surface CPA concentration) is linked to post-thaw CSR contours. Finally, the CSR predictions of this model are compared to those of the previous model that only links the CPA differentials created during loading to postthaw CSR. 2.7 Description of the physical process for the heat transfer model To develop the model to describe the temperature profile during slow cooling and fast thawing, the following assumptions were made. 1. The encapsulated cells are dispersed evenly through the material. 2. The cryoprotecting agents are at a preset concentration and non-permeable solutes at isotonic concentration. 3. There are no temperature differences between a given cell and the immediate extracellular microenvironment during cooling. 4. The concentration of CPA and non-permeable solutes within the construct is in equilibrium with the concentration in the solution outside the beads. 5. A convective heat transfer resistance at the bead surface is neglected and hence the surface temperature is equal to that of the outside solution.

bioartificial tissue construct [10] and adopted here with modification. The heat balance for the control volume is expressed in equation (11) and there is no resistance to heat transfer at the surface of the bead as in equation (12). The environmental temperature cools at a constant rate until equilibrium is achieved with the surroundings (equation (13)). Extracellular space contains many potential ice nucleators and is less likely to supercool than intracellular water (equation (14)). Hence, heterogeneous freezing will predominate (equation (15)) [18]. To predict the rate of ice formation and freeze-concentration, the water–DMSO–NaCl phase diagram is used to predict the liquid water, ice fraction, and solute concentrations [23]. Since experimental data exist for extracellular heterogeneous nucleation, published phase diagrams were used to validate the heterogeneous nucleation rates by linking the solid-fraction–temperature data to the temperature–time data from the heat transfer model [24]. The melting rate (equation (16)) is proportional to the solid phase and is assumed to be linear [10]. Coupled equations for the cell uptake are included for instances where a phase change and cell volume dynamics are incomplete – these equations are the same as equations (5) to (9). Equation (8) is modified to incorporate loss of liquid water via phase change in equation (17). The temperature dependence for the hydraulic conductivity and the solute permeability are modelled as in equations (18) and (19). Lastly, the mass balance for extracellular water fraction is given in equation (20). The model was validated by deploying tissue-specific parameters followed by comparison to experimental data on water transport during freezing of tissue [25] – the settling time and final volume fraction were consistent. (10)

Vi = Vie + Vic rCp

∂T = r  krT + I s L + I shet L ∂t

(12)

2.8 Implementation and validation of the heat transfer model

Ts = Te

The control volume analysis is performed with the control volume consisting of the cell volume (Vc) and the extracellular volume (Ve). The extracellular volume consists of the solution (cryoprotecting agent, salt, and water) and solid phase (biomaterial), while the cell volume is composed of cell water, salts, CPA, and the inactive cell volume. The mathematical models to describe the one dimensional conduction with phase changes have previously been developed for an engineered

Te ðt Þ = Te ðt = 0Þ 

Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

I s = Oexp

(11)



k 2 T3 DTeq

s Ihet = Ohet exp



∂T t and Te ðt = ‘Þ = Tstorage ∂t



khet 2 T3 DTeq

dI  = cDTeq V I dt

(13)

(14) 

(15)

(16)

Predicting the survival rate of mouse embryonic stem cells

   ∂VW Nd VW 0 i =  Lp ART M e   Mn0 ∂t VW VW ∂I s VW ∂t

(17)

   ELp 1 1  Lp = Lpg exp  R T Tref

(18)

   Ea 1 1  Pd ðT Þ = Pdref exp  R T Tref

(19)

nbvw

3

ð

(20)

LINKING THE HEAT TRANSFER MODEL TO FREEZE-CONCENTRATION

Experimental work on freeze-concentration [24] was used to link the temperature model to the freezeconcentration (FC) effect. The temperature–time profile for each radius was modelled and then each mesh-point concentration of the remaining liquid was updated dynamically as the concentration–time profile was generated. Points outside the experimental limit were modelled using theoretical models for freeze-concentration [23]. Finally, the concentration–time curve for the loading, freezeconcentration during cooling, exposure during thawing, and resuspension was integrated to give the exposure in ‘molar-hours’ – this exposure was used to substitute the ‘exposure time’ of the diffusion-membrane transport model in the upcoming section on ‘construction of survival rate contours’. 4

NUMERICAL SOLUTIONS

The mathematical model was implemented in a Matlab (R2008a)TM environment using the PDE solver pdepe. The numerical formula and time step were varied dynamically to suit the nature of the problem. The same integration constraints were used as for the diffusion-membrane transport model. 5

cryopreservation experiments. A cubic smoothing spline algorithm (MatlabÒ) was used to construct survival rate contours based on these data. The smoothing spline has previously been defined as in equations (21) to (23) [26] but is adapted here to the needs of the predictive model. The CSR for each concentration and equilibrium exposure time can be modelled using a non-parametric mean s. For a fixed concentration, CSR depends on exposure time as follows

ð

s = 1  I dt  Ihet dt  Vc=Vi s

CONSTRUCTION OF THE SURVIVAL RATE CONTOURS

Cell survival rate prediction was based on cell survival rate contours drawn on a two-dimensional (2D) domain of equilibrium exposure times and DMSO concentrations. The equilibrium exposure times (alternatively ‘exposure in molar-hours’ for freeze-concentration), DMSO concentrations, and CSR data were derived from suspension

1097

y i = s ðt i Þ + ei

(21)

where ti , i = 1, 2, . . . , n

(22)

The smoothing spline estimator, s, will be the minimizer of equation (23), which minimizes excursions from the experimental data and penalizes roughness n X

ð 00 ½yi  sðti Þ2 + l s ðti Þ2 dt

(23)

i=1

The smoothing spline estimator provides data values that are then plotted in a three-dimensional (3D) contour plot in Matlab. The contour plot interprets each point in the concentration–time domain as a height above the X–Y plane. 6

PREDICTION OF LOCAL AND OVERALL CELL SURVIVAL

The diffusion-membrane transport model provides the concentration and equilibrium exposure time at each point in an alginate bead. The two parameters map on to a particular survival rate contour. This is the per cent local CSR at a given radius; the overall CSR is computed from local cell survival. First, the per cent local CSR was converted to a decimal; then the cell number differential was calculated. The two data were plotted and a curve was fitted to the data. The function describing the curve was then integrated over the domain of the cell number differential yielding the total number of viable cells within the bead. The total number of cells was calculated by multiplying the cell density by the total bead volume. The overall CSR was calculated by dividing the total number of viable cells by the total number of cells in a single bead. The overall CSR was then expressed as a percentage. The overall cell survival rate (OCSR) is mathematically expressed as follows Pn OCSR =

 (r  dBi ) rB

i = 1 yi

(24)

Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

1098

S Sambu, X Xu, H Ye, and Z F Cui

The model implementation is as follows: a curve is fitted on to a viability (y)-cell number differential (r*dB) graph. The resulting function is integrated analytically (Matlab symbolic math toolbox) yielding the integral, which can be evaluated at specified points yielding the total viable cell number. The denominator is calculated, and the resulting quotient (overall CSR) is expressed as a percentage. 7

EXPERIMENTAL MATERIALS AND METHODS

Unless otherwise stated, all materials were sourced from Sigma-Aldrich, UK. Cell culture and preparation: The ESF 183 mESC cell line was derived from mouse strain 129 S2/Sv from the Gardener Lab (University of Oxford, UK) using techniques in the public domain [27]. Mitomycin-treated primary mouse embryonic feeder cells (Department of Zoology, University of Oxford, UK) were initially cultured in Dubelcco’s modified Eagle’s medium (DMEM) supplemented with 10 per cent foetal calf serum (FCS), 1 per cent (v/ v) glutamine, 1 per cent (v/v) beta-mercaptoethanol and 1 per cent (v/v) penicillin/streptomycin. The mESCs were then seeded on to the feeder layer in the mESC culture medium, consisting of DMEM supplemented with 15 per cent (v/v) foetal calf serum supplemented with 1 per cent (v/v) non-essential amino acids, 1 per cent (v/v) b-mercaptoethanol (b-ME), 1 per cent (v/v) penicillin/streptomycin, and mouse leukaemia inhibitory factor (1000 units/ml, Millipore, UK). The mESCs were incubated at 5 per cent vol. CO2 in air and 37 °C. Serum-free and feeder-free stem cell culture: To enrich the stem cells and eliminate feeder cells, the culture method was gradually changed to feederfree culture. Mouse ES cells grown on feeders were passaged and then cultured on a 0.1 per cent gelatine-coated surface in clonal grade medium (CGM, Millipore UK) according to manufacturer’s protocols. Cell encapsulation: Cells were centrifuged then resuspended in 1.2 wt% alginate solution in 0.9 wt% NaCl at a density of 1 3 106 cells/ml. The cell suspension was transferred into a syringe (BD, UK) and ejected, drop by drop, through a 30 gauge needle (BD Bioscience, UK) into a beaker of 102 mM CaCl2 solution. Alginate beads were allowed to polymerize for 10 min before decantation. Excess CaCl2 was removed by washing the alginate beads twice with 0.9 wt% NaCl. The encapsulated cells were then resuspended in culture medium. CPA cryopreservation and thawing: A 50 ml CPA stock solution (20 per cent (v/v) DMSO, 20 per cent Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

(v/v) FBS, 60 per cent (v/v) DMEM) was diluted to the final concentration (10 per cent v/v) in 50 ml of mESC culture medium. The 100 ml, 10 per cent (v/v) DMSO solution contained 10 per cent (v/v) DMSO, 10 per cent (v/v) FBS, 30 per cent (v/v) DMEM, and 50 per cent (v/v) mESC culture medium. DMSO and DMEM concentrations in stock solution were varied when different concentrations of DMSO were required. Cells were cryopreserved in suspension at concentrations of 6, 10, and 12 %vol. DMSO after loading the cryoprotectants for durations ranging from 5 to 60 min at 4 °C, as detailed in Table 2. All samples were cooled at 1 °C/min using a controlledrate freezer (model: Kryo 560-16, Planer plc, UK). Cell suspensions were moved to a 280 ° C freezer after the sample temperature had equilibrated with the chamber temperature at 280 °C. Cryopreserved cell suspensions were rapidly warmed by immersion in a 37 °C water bath until fully thawed. The cryopreservation medium was gradually diluted drop by drop with mESC culture medium immediately after thawing, reducing an original 10 %vol. solution to 2 %vol. The samples were allowed to stand for 10 min at room temperature prior to centrifugation and resuspension in mESC culture medium. For cryopreservation in alginate beads, encapsulated cells were loaded with a 10 %vol. DMSO CPA solution for 20, 30, and 45 min durations at 4 °C. The slow cooling and thawing process for alginate beads is identical to the procedures for suspension cryopreservation. Evaluation of cell survival: To determine cell survival, samples were assayed for cell membrane integrity. Cells embedded in alginate were stained with trypan blue after complete dissolution with filter-sterilized 0.055M sodium citrate in 0.9 %wt NaCl solution (pH 7.4), centrifugation (148g for 5 min), and immediate resuspension in culture medium. Cells in suspension were centrifuged, the supernatant removed and resuspended in 30 uL mESC culture medium. To discern live and dead cells, 30 uL sterile filtered 0.4 %wt/vol. trypan blue was added to the resuspended cells. The number of live cells was quantified using a hemacytometer (Neubauer, UK) within 5 min of staining and expressed as a percentage of total cell number. Statistical analysis: Each experiment was performed in triplicate and repeated three times. Data are presented as mean 6 standard error. Statistical significance between groups was determined by analysis of variance test. Data were presented as means 6 standard error and a p value of less than 0.05 was considered statistically significant.

Predicting the survival rate of mouse embryonic stem cells

Table 2 Experimental parameters and results for suspension cryopreservation. All experiments were conducted using same passage cells cultured under the same conditions. Values represent sample mean and standard error, n = 9 for all data points Time (min)

5 15 20 30 45 60

8

Concentration (%vol.) 6%

8%

10%

12%

3660.6 3660.7 3960.7 4360.5 4060.4 3161.2

4460.4 7160.4 5360.5 4961.2 4560.7 4460.7

4460.6 8160.4 6160.6 5261.2 5060.8 4561.7

3960.6 5060.6 4360.7 3860.9 3861.0 3261.1

RESULTS AND DISCUSSION

8.1 Validation of the diffusion model and membrane transport model The diffusion model was validated by comparison to an analytical solution [28]. The analytical solution is provided in equation (25). The model and analytical solutions are compared using the reference values in Table 1. The results are shown in Fig. 1A. The line graphs for the analytical and numerical solution are coincident showing that there is an agreement between the two techniques. ‘ X

n+1

2j ð1Þ fC ðr, 0Þ  Cs g p n n=1     2 2 1 npr Dn p t sin exp r j j2

C = Cs +

(25)

1099

After validating the diffusion model, the membrane transport model is validated by comparing mouse C57BL/6 ES cell response to CPA loading (1.0M, 295 K) in suspension to the diffusionmembrane transport model. The secondary experimental data for the membrane transport model were obtained from Kashuba Benson et al. [15]. A coefficient of determination (R2) is calculated to determine how well the model fits the data. The result, shown on Fig. 1B, indicates that 91 per cent of the variability in the data is explained by the model. Hence, the model is a good fit for the experimental data. Following this validation, the sensitivity of the model to cell density was examined. Figure 2 shows the impact of cell densities on CPA concentration. From this result, it appears that cell numbers change the CPA concentration transient slightly. The model is improved by taking this aspect into consideration. However, the time (1314 s) taken to achieve the steady-state concentration is virtually unchanged.

8.2

Computations of concentration transients for different model parameters

Following validation and sensitivity analysis, the model was used to derive CPA concentration transients for several parameters. In Fig. 3A, longer loading times lead to higher CPA concentrations. Longer loading times are known to reduce cell viability. This is due to toxicity of CPA during loading [29, 30]

Fig. 1 A two-part validation of the diffusion-membrane transport model. (A) Validation of the numerical diffusion model – the analytical solution (line) and numerical results (triangles) are coincident. (D (4.89 3 10210 m2/s), Hindrance (1.2), Cs (1.4M), R (0.6 mm) and Ti (4 °C)). (B) Validation of the two-parameter (2P) membrane transport model – experimental mESC volume change during CPA loading in suspension (data from Kashuba Benson et al. [15]) are well predicted by the 2P model, R2 = 0.91. Experimental conditions: [DMSO] = 1.0M, hydraulic conductivity, Lp = 0.4160.03um/min/atm, solute permeability, Ps = 4.59 6 0.41 um/min, Vw0 = 50.3 per cent of isosmotic volume (695 um3), and T = 295 K Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

1100

S Sambu, X Xu, H Ye, and Z F Cui

Fig. 2 Simulating the impact of cell density on CPA concentration profile for the bead’s centre (r/ R = 0) and close to the surface (r/R = 0.8). Hindrance varies with cell density (D (4.89 3 10210 m2/s), Cs (1.4M), R (0.6 mm) and Ti (4 °C))

and thawing [31]. For DMSO, this phenomenon varies across cell types and species [29, 30, 32]. After varying the loading time, the impact of the construct radius was explored as depicted in Fig. 3B. Equilibration is faster for smaller construct radii. However, there are operational limitations for constructs at a radius of 0.15 mm or less. They need sophisticated equipment, degrade rapidly in aqueous environments [9], and are susceptible to cryodamage [9, 33]. However, constructs significantly larger than 1.5 mm suffer from larger CPA differentials than constructs below 1.5 mm radius. Furthermore, the effect of varying the hindrance of the CPA in alginate was investigated. Hindrance refers to the impedance to the diffusion of the CPA into alginate. From Fig. 3C, a higher hindrance led to a longer equilibration time because the effective diffusion of the CPA was reduced. In addition to examining the impact of hindrance, the influence of surface concentration (Cs) on the CPA distribution was investigated. From Fig. 3D, increasing the Cs does not alter the time required for the local concentration at the centre of the construct to match the concentration at the surface. While setting the Cs may alter the actual molar concentration of the CPA at various locations within the construct, the time required for equilibration will not change. The heat transfer model prediction on the temperature and concentration transients do show that freezing concentration is significant (see Fig. 4A and Fig. 4C). The remaining concentration of DMSO is increased from 1.4M to over 6M. It has to be pointed out that in the exercise this final concentration is calculated simply based on water–salt–DMSO phase diagram and possible vitrification is completely Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

ignored. In reality when the CPA concentration has reached a certain value, the remaining unfrozen region, excluded by the ice, is simply vitrified. The focus of this paper is not, however, the freeze-concentration, as this is far more complicated to have a thorough study. The purpose is to demonstrate whether freeze-concentration imposes additional heterogeneity to the cells at different location in the beads, hence making the prediction of CSR based on CPA loading parameters invalid. Figure 4 shows the concentration of CPA in the remaining unfrozen liquid during the freezing process (A and B) and thawing (C and D) for two bead sizes. In small constructs (radius at 0.6 mm), the concentration variation induced by freeze-concentration is negligible – i.e. all cells will experience the same concentration transient profile during freezing (see Fig. 4A) and thawing (Fig. 4C). Therefore the main difference would still be the different effect induced during the CPA loading. This is not the case for larger construct (radius 25 mm), the concentration of cells experienced during freezing (Fig. 4B) and thawing (Fig. 4D) depends on the cell location and CPA loading protocol is not the only factor differentiating them any longer. It is hence noted that the analysis developed in the previous sections can only be applied to smaller beads. Fortunately, this would be the case for practical applications.

9

DERIVATION OF THE EQUILIBRIUM EXPOSURE TIME FROM MODEL RESULTS

The DMSO concentration transients have been computed. To connect the model to the experimentally derived survival rate data, the ‘equilibrium exposure time’ will be illustrated using an example. Within the context of 3D cryopreservation, the ‘equilibrium exposure time’ (tE) was defined as the time that cells at a given radius of the alginate bead have been exposed to CPA at a molarity equivalent to Cs. A visual representation of this concept is provided in Fig. 5. To model freeze-concentration, the area under the concentration–time transient is used to measure exposure throughout the loading, cooling, and thawing processes.

10

DERIVATION OF THE SURVIVAL RATE CONTOURS AND THE DETERMINATION OF THE DMSO OPERATION WINDOW

The model can compute the concentration transients for DMSO within the bead and within the cell. To predict CSR, these outputs must be mapped on

Predicting the survival rate of mouse embryonic stem cells

1101

Fig. 3 CPA profiles within an alginate bead under different parameters: (A) loading time, (B) bead radius, (C) hindrances, and (D) surface concentration. For graphs B, C, and D, concentrations are monitored at the centre (r = 0 mm). Other parameters were fixed as per Table 1 (D (4.89 3 10210 m2/s, Ti (4 °C))

to CSR contours based on experimental results for cells cryopreserved in suspension. For purposes of CSR prediction for 3D cryopreservation, a smoothing spline was fitted to the suspension cryopreservation data. Fitting data by splining is useful when the overarching function determining the behaviour of the measured phenomenon is difficult to derive in closed form or is inefficient to handle due to costs in memory and processing speed. Splines guarantee a smoothness of fit through the continuity of the interpolant, its first and second derivatives (including partial derivatives) [34]. Therefore, the ‘spline interpolant’ can be expected to ‘respond’ to the independent variable(s) as the actual ‘overarching function’ would. The smoothing spline model is shown in Fig. 6A (CPA loading only) and 6B (CPA loading and freeze-concentration effects). Part of the experimental objective was to determine the operation window of this cryopreservation system. Formally, the operation window refers to a prescribed time that a population of cells can be exposed to a particular CPA at a given concentration

without a significant loss in viability. Hence, the operation window will have a lower and upper equilibrium exposure time. The operation window is obtained by setting a minimum survival rate contour. For example, by setting the minimum survival rate contour to 60 per cent, the operation window is 12–19 min for 1.4M DMSO (Fig. 6A).

11

COMPARISON OF EXPERIMENTAL DATA FOR SUSPENSION AND 3D CRYOPRESERVATION

Following the derivation of survival rate contours and DMSO operation window, the 3D cryopreservation data are compared to suspension cryopreservation data. The aim is to see the degree of disparity in data. The experimental data for cells cryopreserved in suspension is different from data for cryopreserved encapsulated cells, as expected (Fig. 7). The cell survivals are significantly different in all three loading times. This is because the alginate bead is a diffusion barrier to DMSO. Hence, the concentration transients for encapsulated cells are

Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

1102

S Sambu, X Xu, H Ye, and Z F Cui

7 center mid surface

6

concentration in moles/Liter

concentration in moles/Liter

7

5 4 3 2 1

0

2000

4000 time in seconds

6000

6 5 4 3 2 1

8000

center mid surface

0

20

7 center mid surface

6

concentration in moles/Liter

concentration in moles/Liter

7

5 4 3 2 2000

4000 time in seconds

6000

8000

(B)

center mid surface

6 5 4 3 2 1

0

80

(C)

(A)

1

40 60 time in seconds

0

20

40 60 time in seconds

80

(D)

Fig. 4 CPA concentration at different locations during cooling (A, B) and thawing (C, D) for small (0.6 mm in A and C) and large (2.5 cm in B and D) constructs. (Cs (1.4 M) and Ti (4 °C); other parameters derived from Table 1)

computed transients and matched to the corresponding ‘concentration–time–survival’ contour in the experimental suspension cryopreservation data. This contour specifies the local CSR. 12

Fig. 5 Illustrating the equilibrium exposure time for beads of different radii. The graph depicts normalized concentration transient against the time that constructs have been in the cryopreservation solution. (Hindrance (1.2), D (4.89 3 10210 m2/s), Cs = 1.4M and Ti (4 °C))

different from those of cells in suspension. The diffusion-membrane transport model is able to take this difference into account by predicting the exact concentration transient for encapsulated cells. The equilibrium exposure time can be drawn from the Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

EXPERIMENTALLY USEFUL CORRELATIONS BETWEEN CRYOPRESERVATION PARAMETERS AND LOCAL CELL SURVIVAL

Using the survival rate contour plot, the impact of loading time on local CSR was quantified (Fig. 8A). The low survival rate at the construct surface is associated with the long exposure times to CPA. Similar phenomena have been observed in other cell types, e.g. bovine chondrocytes [31]. Furthermore, the impact of size on local CSR is simulated. The resulting survival rate profile is shown in Fig. 8B. A 1 mm radius construct yielded better survival rates than smaller sized constructs since it had lower equilibrium exposure times. Further, using the concentration profile for different hindrances, the CSR was simulated as shown in Fig. 8C. The survival rate profile improves with increasing hindrance due to

Predicting the survival rate of mouse embryonic stem cells

1103

Fig. 6 Graph A is a 3D, smoothed-spline contour plot based on results from the cryopreserved mESC cell suspensions showing a global maximum at ca. (10 %vol. DMSO, 15 min); graph B shows a similar implementation taking into account the freeze-concentration effects. All experiments were loaded at a temperature of 4 °C, slow cooled at 21 °C/min and stored at 280 °C. Reference values used are provided in Table 1

differentials caused by freeze-concentration. CPA loading differentials are the main reason to explain the differences in post-thaw CSR. 12.1 Comparison of the experimental and predicted cell survival rate

Fig. 7 A comparison of overall cell survival rate for cells cryopreserved in three dimensions and in suspension. All experiments were conducted using 10 %vol. DMSO loaded at 4 °C, slow cooled at 21 °C/min and stored at 280 °C. For 3D cryopreservation, cells were encapsulated in 1.2 wt% alginate

the decrease in equilibrium exposure times. From Fig. 8D, varying Cs away from the optimum (1.4M or 10 %vol. DMSO for mESCs) will lower the CSR within the construct. To compare the impact of incorporating freezeconcentration (FC) into the model, the results in Fig. 8A (not incorporating FC) are compared to the predictions obtained when the exposure to higher concentration of CPA during freeze-concentration is considered (see Fig. 9). From Fig. 9, the two predictions are well matched showing that the FC module does not significantly improve on the mapping between loading and post-thaw CSR provided by the smoothing spline model because small constructs do not suffer significant concentration

After predicting the local CSR, the model was compared to experimental results in 3D cryopreservation. The model and experimental data are compared in Fig. 10. For the 30–45 min loading time, both results show a decreasing overall CSR with increase in loading time. However, for the 20–30 min loading time, the trends are in conflict – while the model predicts a slight decrease in overall cell survival, the experiments showed a sharp rise. Coupled with the dissimilarity in absolute values for the overall cell survival, these results show that there are additional factors affecting the overall CSRs. Such factors may include the nature of the encapsulation material [3], the range of bead sizes, and biases in cell location within the construct. A bias in cell location causes a varying CSR across the construct (in the simulations) and can be observed in situ via 3D fluorescence or multiphoton microscopy (data not provided). Nonetheless, for the best performing loading time (30 min), the predicted and experimental results are in agreement. 13

PREDICTION OF THE OVERALL CELL SURVIVAL RATE

After confirming that the diffusion-membrane transport model could predict CSR in 3D cryopreservation, it was used to simulate the effect of changing a range of parameters on cell survival. The Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

1104

S Sambu, X Xu, H Ye, and Z F Cui

Fig. 8 The predicted local cell survival rate profile for different model parameters based on concentration–time profiles and survival rate contour plot: (A) loading time, (B) bead radius, (C) hindrances, and (D) surface concentration

hindrances can be caused by high cell densities [35]. This will reduce equilibrium exposure times and decrease cell viability. Hence, a trade-off needs to be made between these competing aims. Lastly, the effect of surface CPA concentrations is shown in Fig. 11D. For DMSO, concentrations lower than 1.4M will lead to losses in viability. This is consistent with previously reported results for other cell types [36]. Concentrations greater than 1.4M also have diminished post-thaw CSR. Overall, suboptimal surface concentrations are unattractive as they lower local and overall cell viability. Fig. 9 A comparison of the diffusion model before the incorporation of the freeze-concentration effect at various loading times and after the incorporation of the freeze-concentration effect (FC). Model parameters: D (4.89 3 10210 m2/s) and Ti (4 °C); other reference parameters were fixed as per Table 1

local CSR data were used to compute the overall CSR for each set of parameters as described in the methods section. For loading time, the overall per cent CSR for the entire construct shows that there exists an operation window within which yields are higher than at all other times – Fig. 11A. The construct radius also affects overall CSR similarly: Figure 11B shows that the overall cell survival (OCS) rate peaks at a 1 mm radius then begins to drop. Furthermore, the overall CSR peaks at a hindrance of 3 then drops thereafter in a fashion similar to loading time and construct size (Fig. 11C). One reason for this drop in the OCS may be the increase in ice fracturing [9]. High construct Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

14

CONCLUSIONS

A numerical model has been implemented to simulate the mass transfer and membrane transport of CPA into a cell-seeded spherical construct. This model can predict the local concentrations,

Fig. 10 Comparison of the model predictions for overall cell survival rate to the experimental results. (Hindrance (1.2), R (0.6 mm), D (4.89 3 10210 m2/s), Cs = 1.4M, and Ti (4 ° C))

Predicting the survival rate of mouse embryonic stem cells

1105

Fig. 11 The predicted overall cell survival rate profile for different model parameters based on concentration profiles and survival rate contour plot: (A) loading time, (B) bead radius, (C) hindrances, and (D) surface concentration. Model parameters: Hindrance (1.2), R (0.6 mm), D (4.89 3 10210 m2/s), Cs = 1.4M and Ti (4 °C). Other parameters were fixed as per Table 1

intracellular CPA concentrations and equilibrium exposure times at different points within the construct. Predictions of concentration transients can be made for different cryopreservation parameters such as loading time, construct radius, hindrances, and surface concentrations. Using experimental suspension data, the relationship between equilibrium exposure time, DMSO concentration and CSR was depicted using a smoothing spline. The survival rate contour plot can be used to derive an operation window for a given CPA. For DMSO, the results indicated an operation window with upper and lower equilibrium exposure times of 12–19 min for a minimum 60 per cent CSR. The survival rate contour plot was then used to predict local and overall CSR based on model results. The model can be applied to various encapsulation materials, CPAs and multistep loading processes. Since these predictions are well matched to predictions from a model incorporating freeze-concentration dynamics for small constructs, the diffusion-membrane transport model is sufficient to predict expected post-thaw CSR for small constructs. The predictions of CSR are closest to the experimental data in the 3D cryopreservation experiments with the highest cell survival. This investigation shows that cell survival-based optimisation is possible using mathematical modelling and experimentally determined parameters.

FUNDING Part of this work was supported by the Biotechnology and Biological Science Research Council, UK.

ACKNOWLEDGEMENTS Sammy Sambu is grateful to the Rhodes Trust for providing the Rhodes scholarship. Ó Authors 2011

REFERENCES 1 Wu, Y., Yu, H., Chang, S., Magalhaes, R., and Kuleshova, L. L. Vitreous cryopreservation of cell– biomaterial constructs involving encapsulated hepatocytes. Tissue Engng, 2007, 13, 649–658. 2 Kampf, N. The use of polymers for coating of cells. Polym. Adv. Technol., 2002, 13, 896–905. ¨ller, K. J., Sukhorukov, V. L., 3 Shirakashi, R., Mu and Zimmermann, U. Effects of physiological isotonic cryoprotectants on living cells during the freezing–thawing process and effects of their uptake by electroporation: Sp 2 cells in alginate– trehalose solutions. Heat Transfer Asian Res., 2003, 32, 511–523. 4 Siti-Ismail, N., Bishop, A. E., Polak, J. M., and Mantalaris, A. The benefit of human embryonic stem cell encapsulation for prolonged feeder-free maintenance. Biomaterials, 2008, 29, 3946–3952. 5 Lim, S. Microencapsulated islets as bioartificial endocrine pancreas. Science, 1980, 210, 908–910. 6 Kopstad, G. and Elgsaeter, A. Theoretical analysis of the ice crystal size distribution in frozen aqueous specimens. Biophys. J., 1982, 40, 155–161. 7 Bhakta, G., Lee, K. H., Magalha˜es, R., Wen, F., Gouk, S. S., Hutmacher, D. W., and Kuleshova, L. L. Cryoreservation of alginate–fibrin beads involving bone marrow derived mesenchymal stromal cells by vitrification. Biomaterials, 2009, 30, 336– 343. Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

1106

S Sambu, X Xu, H Ye, and Z F Cui

8 Fahy, G. M., MacFarlane, D. R., Angell, C. A., and Meryman, H. T. Vitrification as an approach to cryopreservation. Cryobiology, 1984, 21, 407–426. 9 Heng, B. C., Yu, Y. -H., and Ng, S. C. Slow-cooling protocols for microcapsule cryopreservation. J. Microencapsulation, 2004, 21, 455–467. 10 Cui, Z. F., Dykhuizen, R. C., Nerem, R. M., and Sembanis, A. Modeling of cryopreservation of engineered tissues with one-dimensional geometry. Biotechnol. Prog., 2002, 18, 354–361. 11 Devireddy, R. V. Predicted permeability parameters of human ovarian tissue cells to various cryoprotectants and water. Molecular Reproductive Dev., 2005, 70, 333–343. 12 Doran, S. J. and Decorps, M. A robust, single-shot method for measuring diffusion coefficients using the ‘‘burst’’ sequence. J. Magn. Resonance A, 1995, 117, 311–316. 13 Green, D. W., Maloney, J. O., and Perry, R. H. Perry’s chemical engineers’ handbook, 1984 (McGrawHill, London). 14 Axelsson, A. and Persson, B. Determination of effective diffusion coefficients in calcium alginate gel plates with varying yeast cell content. Appl. Biochem. Biotechnol., 1988, 18, 231–250. 15 Kashuba Benson, C. M., Benson, J. D., and Critser, J. K. An improved cryopreservation method for a mouse embryonic stem cell line. Cryobiology, 2008, 56, 120–130. 16 Blandino, A., Macias, M., and Cantero, D. Formation of calcium alginate gel capsules: Influence of sodium alginate and CaCl2 concentration on gelation kinetics. J. Biosci. Bioengng, 1999, 88, 686–689. 17 Yarmush, M. L., Toner, M., Dunn, J. C. Y., Rotem, A., Hubel, A., and Tompkins, R. G. Hepatic tissue engineering: Development of critical technologies. Ann. N. Y. Acad. Sci., 1992, 665, 238–252. 18 Toner, M., Cravalho, E. G., and Karel, M. Thermodynamics and kinetics of intracellular ice formation during freezing of biological cells. J. Appl. Phys., 1990, 67, 1582–1593. 19 Kleinhans, F. W. Membrane permeability modeling: Kedem–Katchalsky vs a two-parameter formalism. Cryobiology, 1998, 37, 271–289. 20 Skeel, R. D. and Berzins, M. A method for the spatial discretization of parabolic equations in one space variable. SIAM J. Sci. Stat. Comput., 1990, 11, 1–32. 21 Shampine, L. F. and Reichelt, M. W. The MATLAB ode suite. SIAM J. Sci. Stat. Comput., 1997, 18, 1– 22. 22 Shampine, L. F., Reichelt, M. W., and Kierzenka, J. A. Solving index-I DAEs in MATLAB and Simulink. SIAM Rev., 1999, 41, 538–552. 23 Pegg, D. E. Equations for obtaining melting-points and eutectic temperatures for the ternary-system dimethylsulfoxide sodium-chloride water. Cryo. Lett., 1986, 7, 387–394. 24 Cocks, F. H. and Brower, W. E. Phase-diagram relationships in cryobiology. Cryobiology, 1974, 11, 340–358.

Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine

25 Han, B. and Bischof, J. C. Engineering challenges in tissue preservation. In Vitro Cell Dev. Biol. Anim., 2004, 40, 11A–11A. 26 Wu, H. and Zhang, J. T. Nonparametric regression methods for longitudinal data analysis, 2006, pp. 149–150 (Wiley-Blackwell, New Jersey, USA). 27 Brook, F. A. and Gardner, R. L. The origin and efficient derivation of embryonic stem cells in the mouse. Proc. Natl Acad. Sci. USA, 1997, 94, 5709– 5712. 28 Perry, J. H., Crank, J., Newman, A. B., Vinograde, B., Maple, C. G., and Holl, D. L. Unsteady state molecular diffusion. In Diffusional mass transfer (Skelland A. H. P., ed.), 1974, pp. 22– 48 (Wiley, New York). 29 Da Violante, G., Zerrouk, N., Richard, I., Provot, G., Chaumeil, J. C., and Arnaud, P. Evaluation of the cytotoxicity effect of dimethyl sulfoxide (DMSO) on Caco2/TC7 colon tumor cell cultures. Biol. Pharm. Bull., 2002, 25, 1600–1603. 30 Qi, W., Ding, D., and Salvi, R. J. Cytotoxic effects of dimethyl sulphoxide (DMSO) on cochlear organotypic cultures. Hear. Res., 2008, 236, 52–60. 31 Tomford, W. W., Fredericks, G. R., and Mankin, H. J. Studies on cryopreservation of articular cartilage chondrocytes. J. Bone Joint Surg., 1984, 66, 253–259. 32 Wusteman, M. C., Armitage, W. J., Wang, L., Busza, A. L., and Pegg, D. E. Cryopreservation studies with porcine corneas. Curr. Eye Res., 1999, 19, 228–233. 33 Chia, S., Leong, K. W., Li, J., Xu, X., Zeng, K., Er, P., Gao, S., and Yu, H. Hepatocyte encapsulation for enhanced cellular functions. Tissue Engng, 2000, 6, 481–495. 34 De Boor, C. A practical guide to splines, revised edition, 2001 (Springer, New York). 35 Pu, H. T. and Yang, R. Y. K. Diffusion of sucrose and yohimbine in calcium alginate gel beads with or without entrapped plants cells. Biotechnol. Bioengng, 1988, 32, 891–896. 36 He, S. and Woods, III, L. C. Effects of dimethyl sulfoxide and glycine on cryopreservation induced damage of plasma membranes and mitochondria to striped bass (Morone saxatilis) sperm. Cryobiology, 2004, 48, 254–262.

APPENDIX Notation A B c C Cp Cs D

cell area (m2) total volume of one alginate bead (m3) melting coefficient CPA concentration (moles/l) heat capacity (J/kg K) CPA concentration at the surface of the alginate bead during CPA loading (moles/l) diffusivity (m2/s)

Predicting the survival rate of mouse embryonic stem cells

ei , ej

ELP Ea H k Lp Me d Me Min0 Nd Pd r R

disturbance (difference between a sample measurement and a non-parametric population mean – ei and ej have mean 0 and constant variance activation energy for hydraulic conductivity (J/mole) activation energy for solute permeability (DMSO) (J/mole) hindrance thermal conductivity (W/m K) hydraulic conductivity (m/Pa per s) osmolality of extracellular CPA (Osmoles/ kg) total osmolality of extracellular solutes (Osmoles/kg) initial osmolality of intracellular nonpermeable solutes (Osmoles/kg) moles of CPA (moles) CPA permeability (m/s) radial distance from construct centre (m) gas constant (J/m per K)

s t T Ti Vb Vc Vd Vd Vi VW0 VW V0 y j l r rm y

1107

smoothing spline estimator time (s) temperature (K) initial temperature (K) volume of cell solids and non-osmotic water (m3) total cell volume (m3) DMSO volume (m3) partial molar volume of CPA (m3/mole) instantaneous volume (m3) initial cell water volume (m3) cell water volume (m3) initial volume (m3) cell survival rate construct radius (m) penalty for roughness (varies between 0 and 1) cell density (cells/m3) material density (kg/m3) local cell survival rate

Proc. IMechE Vol. 225 Part H: J. Engineering in Medicine