Computational Modeling for in vitro Tissue Cultivation Kyriacos - - PDF document

Computational Modeling for in vitro Tissue Cultivation

Kyriacos Zygourakis

Chemical and Biomolecular Engineering Rice University

August 12, 2015

Scales and Modeling Approaches for Biological Systems

Range of time scales is even broader (10-9-109 s)

Models for in vitro Tissue Cultivation

Scaffold Design (Biomimetics)

? ? ?

Single Cell Models Bioreactor Design Growth Factors

Cells migrate, proliferate and differentiate

Tissue Growth Modulated by Mass Transfer Processes

• Decreased availability of 


nutrients (and growth factors) in 
 the interior of 3D scaffolds at
 high cell densities.

• The viable size of bioartificial 


constructs is limited to a few 
 hundred microns due to hypoxia, 
 nutrient insufficiency and/or 
 waste accumulation

• Large tissue replacements fail due to


necrosis at the central region.


(Sikavitsas et al., J. Biomed Mater. Res. ,

62: 136-148, 2002)

Strategies to enhance vascularization

• Directing cell behavior through

growth factor delivery

• Using co-culturing systems
• Engineering biomaterials with

appropriate properties (biomimetics)

• Incorporating microfabrication

techniques

• Applying mechanical

stimulation (when necessary)
 


(Khademhosseini et al, Tissue Eng, 2012)

Initial and Boundary Conditions Are Important!

Scaffold Design (Biomimetics)

? ? ?

Single Cell Models Bioreactor Design Growth Factors Seed cell distribution

Signals

Competition of Dynamic Processes

Cells

Scaffolds Dynamics Cell Population Dynamics

• Migration
• Proliferation
• Differentiation

Mass Transfer Diffusion and Uptake of Nutrients and Growth Factors ICs & BCs

• Bioreactor
• Scaffold
• Seed cell distribution

Culture Media Culture Media

z = 0 z = L

Nutrient Diffusion and Consumption

∂C ∂t = De ∂2C ∂z2 − ρcell vmaxC Km + C

Diffusion coefficient Diffusion Term Cell density Nutrient consumption rate

Cs Cs

Nutrient Diffusion and Consumption

Base case parameters: Temporal evolution of nutrient concentration profiles ∂C ∂t = De ∂2C ∂z2 − ρcell vmaxC Km + C C 0

( ) = Cs

C L

( ) = Cs

L = 0.001 m

De = 7 ×10−11 m2/s ρcell ⋅vmax = 8.3×10−4 mol/(m3 ⋅s) Cs = 5 mM

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Nutrient concentration, mol/m3 Normalized Distance, z/L Initial 45 min 90 min 135 min 180 min mM

Nutrient Diffusion and Consumption

Base case parameters: ∂C ∂t = De ∂2C ∂z2 − ρcell vmaxC Km + C C 0

( ) = Cs

C L

( ) = Cs

L = 0.002 m

De = 7 ×10−11 m2/s ρcell ⋅vmax = 8.3×10−4 mol/(m3 ⋅s) Cs = 5 mM

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Nutrient concentration, mol/m3 Normalized Distance, z/L Initial 60 min 120 min 180 min 240 min

Temporal evolution of nutrient concentration profiles

mM

Nutrient Diffusion and Consumption

∂C ∂t = De ∂2C ∂z2 − ρcell vmaxC Km + C C 0

( ) = Cs

C L

( ) = Cs

L = 0.004 m

De = 7 ×10−11 m2/s ρcell ⋅vmax = 8.3×10−4 mol/(m3 ⋅s) Cs = 5 mM

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Nutrient concentration, mol/m3 Normalized Distance, z/L Initial 60 min 120 min 180 min 240 min

Temporal evolution of nutrient concentration profiles Base case parameters:

mM

Effect of maximum nutrient uptake rate

∂C ∂t = De ∂2C ∂z2 − ρcell vmaxC Km + C C 0

( ) = Cs

C L

( ) = Cs

Base case parameters: Diffusional limitations becomes stronger when the uptake rate constant vmax increases

L = 0.002 m De = 7 ×10−11 m2/s ρcell ⋅vmax = 8.3×10−4 mol/(m3 ⋅s) Cs = 5 mM

1 2 3 4 5 0.25 0.5 0.75 1 1.25 1.5 1.75 2

Nutrient concentration, mM Distance, mm 2 vmax vmax 0.5 vmax

1 2 3 4 5 0.25 0.5 0.75 1 1.25 1.5 1.75 2

Nutrient concentration, mM Distance, mm 10 De De 0.1 De

Effect of nutrient diffusion coefficient

Base case parameters: Diffusional limitations depend strongly on the magnitude of effective diffusivities

L = 0.002 m De = 7 ×10−11 m2/s ρcell ⋅vmax = 8.3×10−4 mol/(m3 ⋅s) Cs = 5 mM Effective diffusivity

∂C ∂t = De ∂2C ∂z2 − ρcell vmaxC Km + C C 0

( ) = Cs

C L

( ) = Cs

1 2 3 4 5 6 7 8 9 10 0.5 1 1.5 2 2.5 3 3.5 4

Nutrient concentration, mM Distance, mm 10 mM 5 mM 2.5 mM

Effect of surface nutrient concentration

Other parameters: High surface concentrations may raise intra- tissue concentrations above the desired levels

L = 0.004 m De = 7 ×10−11 m2/s ρcell ⋅vmax = 8.3×10−4 mol/(m3 ⋅s) Surface nutrient concentration

∂C ∂t = De ∂2C ∂z2 − ρcell vmaxC Km + C C 0

( ) = Cs

C L

( ) = Cs

How to quickly estimate the extent of diffusional limitations...

The Thiele Modulus

Introduce dimensionless variables: u = C Cs , ζ = z L , τ = Det L2 to obtain: ∂u ∂τ = ∂2u ∂ζ 2 − φ 2 u β + u 0<ζ <1, τ >0 u 0,τ

( ) = u 1,τ ( ) = 1 and u ζ,0 ( ) = u0 ζ ( )

φ = L ρcell ⋅vmax De ⋅Cs Thiele modulus= Consumption rate Diffusion rate β= Km Cs

The Thiele Modulus

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Normalized nutrient concentration, C/Cs Normalized Distance, z/L

φ = 0.77 φ = 1.54 φ = 3.08 φ = 6.15 φ = 12.3 ϕ ↑ ⇓ penetration depth ↓

3D Models of Ewing Sarcoma Tumors

SEM micrograph of human EWS cells seeded in electrospun 3D PCL scaffold Fong et al, PNAS (2013)

Response of human EWS cells to doxorubicin φ = L ρcell ⋅vmax De ⋅Cs Cs ↑ ⇒ ϕ ↓ ⇒ penetration depth ↑

3D Models of Ewing Sarcoma Tumors

SEM micrograph of human EWS cells seeded in electrospun 3D PCL scaffold Fong et al, PNAS (2013)

Response of human EWS cells to doxorubicin φ = L ρcell ⋅vmax De ⋅Cs vmax ↓ ⇒ ϕ ↓ ⇒ penetration depth ↑

3D Models of Ewing Sarcoma Tumors

φ = L ρcell ⋅vmax De ⋅Cs ρcell ↑ ⇒ ϕ ↑ ⇒ penetration depth ↓

SEM micrograph of human EWS cells seeded in electrospun 3D PCL scaffold Fong et al, PNAS (2013)

Response of human EWS cells to doxorubicin

How can we overcome diffusional limitations?

Perfusion Bioreactor Systems

Bancroft, Sikavitsas and Mikos, Tissue Engineering, 9, 549-554 (2003)

Continuous flow of media through the scaffold

Convection, Dispersion and Consumption

∂C ∂t + vz ∂C ∂z = De,z ∂2C ∂z2 − ρcellVmaxC Km + C

Convection Consumption Dispersion

Boundary conditions: De,z ∂C ∂z = vz C0 − C

( ) at z = 0

∂C ∂z = 0 at z = L

z = L z = 0

Assumptions:

• Axial (1-D) flow
• No concentration variations in 


radial direction

• Uniform cell concentration

C0

Convection, Dispersion and Consumption

∂u ∂τ + Pe ∂u ∂ζ = ∂2u ∂ζ 2 − φ 2 u β + u

Dimensionless variables: u = C Cb , ζ = z L , τ = Dzt L2 Dimensionless numbers: φ 2 = L2 ρcell ⋅vmax Dz ⋅Cb , Pe = L ⋅Vz Dz , β= Km Cb

Peclet Number Convection Thiele modulus (Damkoeler number) Consumption Dispersion

Boundary conditions: ∂u ∂ζ = Pe u −1

( )

at ζ = 0 ∂u ∂ζ = 0 at ζ = 1

Will perfusion overcome the diffusional limitations?

Base case parameters: ∂C ∂t = De ∂2C ∂z2 − ρcell vmaxC Km + C C 0

( ) = Cs

C L

( ) = Cs

De = 7 ×10−11 m2/s ρcell ⋅vmax = 8.3×10−4 mol/(m3 ⋅s) Cs = 5 mM

Pure Diffusion Steady-state concentration profile

Normalized Distance, z/L

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Nutrient Concentration, mM

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

L = 2 mm

Performance of perfusion bioreactors

∂u ∂τ + Pe ∂u ∂ζ = ∂2u ∂ζ 2 − φ 2 u β + u ∂u ∂ζ = Pe u −1

( )

at ζ = 0 ∂u ∂ζ = 0 at ζ = 1 vz = 6.2 ×10−6 m/s ρcell ⋅vmax = 8.3×10−4 mol/(m3 ⋅s) C0 = 5 mM

Normalized Distance, z/L

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Nutrient Concentration, mM

1 2 3 4 5 6

1 min 2 min 3 min 4 min 5 min 6 min Steady-state profile

Base case parameters: Perfusion can help maintain nutrient concentration constant across the tissue construct

L = 2 mm

4.8 mM

Normalized Distance, z/L

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Nutrient Concentration, mM

1 2 3 4 5 6

Performance of perfusion bioreactors

∂u ∂τ + Pe ∂u ∂ζ = ∂2u ∂ζ 2 − φ 2 u β + u ∂u ∂ζ = Pe u −1

( )

at ζ = 0 ∂u ∂ζ = 0 at ζ = 1

2 min 4 min 6 min 8 min 10 min 12 min Steady-state profile

Perfusion can help maintain nutrient concentration constant across the tissue construct

vz = 6.2 ×10−6 m/s ρcell ⋅vmax = 8.3×10−4 mol/(m3 ⋅s) C0 = 5 mM

Base case parameters:

L = 4 mm

4.6 mM

Normalized Distance, z/L

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Nutrient Concentration, mM

1 2 3 4 5 6

Performance of perfusion bioreactors

∂u ∂τ + Pe ∂u ∂ζ = ∂2u ∂ζ 2 − φ 2 u β + u ∂u ∂ζ = Pe u −1

( )

at ζ = 0 ∂u ∂ζ = 0 at ζ = 1

4 min 8 min 12 min 16 min 20 min Steady-state profile

Perfusion can help maintain nutrient concentration fairly constant across the tissue construct

vz = 6.2 ×10−6 m/s ρcell ⋅vmax = 8.3×10−4 mol/(m3 ⋅s) C0 = 5 mM

Base case parameters:

L = 8 mm

4.3 mM

Performance of perfusion bioreactors

High cell densities or high nutrient consumption rates decrease the effectiveness of perfusion. Nutrient transport rates cannot keep up with consumption!

vz* = 6.2 ×10−6 m/s ρcell ⋅vmax

( )* = 8.3×10−4 mol/(m3 ⋅s)

C0* = 5 mM

Normalized Distance, z/L

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Nutrient Concentration, mM

1 2 3 4 5 6

ρcell ⋅Vmax

( )*

20 × ρcell ⋅Vmax

( )*

50 × ρcell ⋅Vmax

( )*

200 × ρcell ⋅Vmax

( )*

L = 2 mm

Parameters:

Normalized Distance, z/L

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Nutrient Concentration, mM

1 2 3 4 5 6

Performance of perfusion bioreactors

Increasing the flow rate improves the effectiveness of perfusion, but exposes the cells to higher shear stresses!

vz* = 6.2 ×10−6 m/s ρcell ⋅vmax

( ) = 1.7 ×10−1 mol/(m3 ⋅s)

C0* = 5 mM

Parameters:

L = 2 mm

vz * 3× vz * 10 × vz * 30 × vz *

How can we reduce the shear stress?

Scaffolds with Perfusion Channels

z = L z = 0

• Forced convection through the channels
• Diffusion of nutrient from the channel walls into the tissue

Scaffolds with Perfusion Channels

Representative Unit Element L d

0.25 0.5 0.75 1 0.25 0.5 0.75 1

x y Finite Element Mesh

∂u ∂τ = ∂2u ∂x2 + ∂2u ∂y2 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ − φ 2 u β + u u = 1

• n the channel wall

∂u ∂x = 0 or ∂u ∂y = 0 on the edges

Scaffolds with Perfusion Channels

0.5 1 0.5 1

x y

τ =1000 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

Plot of C x,y,t

( ) = u ⋅Cs

φ 2 = 3

∂u ∂τ = ∂2u ∂ξ 2 + ∂2u ∂η2 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ − φ 2u u = 1

• n ∂Ωchannel

∂u ∂n = 0

• n ∂Ωedges

Scaffolds with Perfusion Channels

φ 2 = 3

Plot of C x,y,t

( ) = u ⋅Cs

∂u ∂τ = ∂2u ∂ξ 2 + ∂2u ∂η2 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ − φ 2u u = 1

• n ∂Ωchannel

∂u ∂n = 0

• n ∂Ωedges

Scaffolds with Perfusion Channels

φ 2 = 30

Plot of C x,y,t

( ) = u ⋅Cs

∂u ∂τ = ∂2u ∂ξ 2 + ∂2u ∂η2 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ − φ 2u u = 1

• n ∂Ωchannel

∂u ∂n = 0

• n ∂Ωedges

But, can we treat the cells as a continuum?

G1 S G2 M

Discrete model for cell population dynamics

Cells are seeded in 3-D scaffolds that allow them to migrate in all directions. Model Parameters: 
 Migration speed S, Persistence P, Division time td, Initial seeding density c0, Directional probabilities pi, i = 1,2,3,4,5,6

Cheng et al., Biophysical J. (2006)

Cells execute persistent random walks, collide and proliferate ...

Discrete model for cell population dynamics

Cells execute persistent random walks, collide and proliferate ... … until we reach confluence.

0.2 0.4 0.6 0.8 1 2 4 6 8 10 12

Cell fraction Time, days

No contact inhibition

Simulation

S = 10 µm/hr, P = 1 hr, c0 = 0.1% Cells are seeded in 3-D scaffolds that allow them to migrate in all directions. Model Parameters: 
 Migration speed S, Persistence P, Division time td, Initial seeding density c0, Directional probabilities pi, i = 1,2,3,4,5,6

Cheng et al., Biophysical J. (2006)

Modulation of Cell Functions

1.5 1.6 1.7 1.8 1.9 2.0 2.1 0.0 2.0 4.0 6.0 Doubling rate r

g

Glucose Concentration, mg/ml

r

g = kgCglu

K + Cglu

The doubling time rg of human diploid fibroblast (HDF) cells shows a Monod-like dependence

• n the glucose concentration in

the culture media.

Cheng and Zygourakis, 2005

Proliferation

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 6 12 18 24 Cell Migration Speed, µm/min Time, hr

Serum-containing media Serum-free media

Changes in the migration speed

• f rat prostate cancer cells

cultured in serum-containing and serum-free media.

Kouvroukoglou et al., 1998

Migration

Hybrid, Multi-Scale Models

10-7 10-6 10-5 10-4 10-3 10-2 10-1 100 Length Scale (m)

Diffusion-Reaction Equation(s): Transport and consumption of nutrients and growth factors in cellularized scaffold Heterogeneous Cell Population Dynamics: Migration, proliferation, cell-cell interactions and cell differentiation Single Cell Models: Intracellular processes modulating cell 
 function (proliferation,migration) Cells Scaffold Bioreactor

Algebraic Equations or Systems of ODEs Discrete Model (CA) BCs Difgusion- Reaction PDEs

Model Components Length Scale (m)

10-7 10-6 10-5 10-4 10-3 10-2 10-1 100

• Transport and consumption of nutrients or growth factors:

The Diffusion-Reaction Problem

∂c ∂t = ∂ ∂x De ∂c ∂x ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ + ∂ ∂y De ∂c ∂y ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ + ∂ ∂z De ∂c ∂z ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ − ρ x,y,z

( )R c ( ) in Ω

c = cs

• n ∂Ω Dirichlet

∂c ∂n = 0

• n ∂Ω Neumann

kg cb − cs

( ) = De

∂c ∂n on ∂Ω Mixed Cheng et al., Biophysical J. (2009)

The Diffusion-Reaction Problem

• Nutrient diffusion and uptake:
• Dimensionless form:

∂c ∂t = ∂ ∂x De ∂c ∂x ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ + ∂ ∂y De ∂c ∂y ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ + ∂ ∂z De ∂c ∂z ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ − ρ x,y,z

( ) Vmaxc

Km + c in Ω ∂u ∂τ = ∂ ∂ξ δ ∂u ∂ξ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ + ∂ ∂ψ δ ∂c ∂ψ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ + ∂ ∂ζ δ ∂c ∂ζ ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ − ϕ 2u β + u in Ω u = c c* , τ = t ⋅ De

*

L2 , ξ = x L , ψ = y L , ζ = z L , δ = De De

* , β = Km

c* Thiele Modulus: φ 2 = L2 ρcellVmax De

*c*

Magnitude of Transport Limitations

Thiele modulus: φ = L ρcellVmax De,tCb Biot Number: Bi = kgL De,t

Monitor cell population and nutrient concentration on

• n central cross section

External mass transfer Internal mass transfer (diffusion) ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ ⎟

Effects of mass transport limitations

Initial

t=0 t=2.5 days t=5 days φ=1.15

Cells

Nutrient Concentration

φ=11.5

Cells

Nutrient Concentration

Effects of mass transport limitations

Initial

t=0 t=2.5 days t=5 days φ=1.15

Cells

Nutrient Concentration

φ=115

Cells

Nutrient Concentration

Experiment and model predictions

(Sikavitsas et al., J. Biomed Mater. Res. ,

62: 136-148, 2002)

Model Predictions

2 mm Cross section of cubic scaffold

Experiment

Heterogeneity and Transport Limitations

C = Cs on ∂Ω (all faces of cube)

• Two sub-populations
• High metabolic rates
• Strong transport

limitations

• Various initial cell

distributions (mixed and segregated) Boundary conditions:

Transport Limitations and Cell Death

Population 1: κ0=0.005, S=50 μm/hr, td=17.4 hr, ρcVmax=1.15x10-3 (mole/(m3)(s), c*/c0=0.008 Population 2: κ0=0.005, S=5 μm/hr, td=22.8 hr, ρcVmax=5.75x10-4 (mole/(m3)(s), c*/c0=0.004

0.00 0.20 0.40 0.60 0.80 1.00 7 14 21

Total Cell Volume Fraction κ(t) Time, days

No Transport Limitations

Mass Transport Limitations No Cell Death With Cell Death

The dynamics of competing biological and physical processes modulate both the structure and the growth rates of tissues.