This tutorial is automatically generated from the file test/python/tutorials//TestPythonLatticeBasedAngiogenesisTutorial.py.

In [1]:
# Jupyter notebook specific imports 
import matplotlib as mpl 
from IPython import display 
%matplotlib inline

A Lattice Based Angiogenesis Tutorial

This tutorial introduces a lattice based angiogenesis problem based on a simplified version of the vascular tumour application described in Owen et al. 2011.

It is a 2D simulation using cellular automaton for cells, a regular grid for vessel movement and the same grid for the solution of partial differential equations for oxygen and VEGF transport using the finite difference method.

The Test

In [2]:
import chaste # Core Chaste functionality
import chaste.cell_based # Chaste Cell Populations
chaste.init() # Initialize MPI and PETSc
import microvessel_chaste # Core Microvessel Chaste functionality
import microvessel_chaste.geometry # Geometry tools
import microvessel_chaste.mesh # Meshing
import microvessel_chaste.population.vessel # Vessel tools
import microvessel_chaste.pde # PDE and solvers
import microvessel_chaste.simulation # Flow and angiogenesis solvers
import microvessel_chaste.visualization # Visualization
from microvessel_chaste.utility import * # Dimensional analysis: bring in all units for convenience
# Set up the test 
chaste.cell_based.SetupNotebookTest()

Set up output file management and seed the random number generator.

In [3]:
file_handler = chaste.core.OutputFileHandler("Python/TestLatticeBasedAngiogenesisTutorial")
chaste.core.RandomNumberGenerator.Instance().Reseed(12345)

This component uses explicit dimensions for all quantities, but interfaces with solvers which take non-dimensional inputs. The BaseUnits singleton takes time, length and mass reference scales to allow non-dimensionalisation when sending quantities to external solvers and re-dimensionalisation of results. For our purposes microns for length and hours for time are suitable base units.

In [4]:
reference_length = 1.e-6 * metre()
reference_time = 3600.0 * second()
BaseUnits.Instance().SetReferenceLengthScale(reference_length)
BaseUnits.Instance().SetReferenceTimeScale(reference_time)

Set up the lattice (grid), we will use the same dimensions as Owen et al. 2011. Note that we are using hard-coded parameters from that paper. You can see the values by printing. When we get the value of the parameter by doing Owen11Parameters.mpLatticeSpacing.GetValue("User") a record is kept that this parameter has been used in the simulation. A record of all parameters used in a simulation can be dumped to file on completion, as will be shown below.

In [5]:
grid = microvessel_chaste.mesh.RegularGrid2()
grid_spacing = Owen11Parameters.mpLatticeSpacing.GetValue("User")
grid.SetSpacing(grid_spacing)
grid.SetExtents([51, 51, 1])

We can write and visualize the grid.

In [6]:
grid.Write(file_handler)
scene = microvessel_chaste.visualization.MicrovesselVtkScene2()
scene.SetRegularGrid(grid)
scene.GetRegularGridActorGenerator().SetVolumeOpacity(0.1)
nb_manager = microvessel_chaste.visualization.JupyterNotebookManager()
nb_manager.vtk_show(scene, height=600, width = 1000)
Out[6]:

Next, set up the vessel network, this will initially consist of two, large counter-flowing vessels. Also set the inlet and outlet pressures and flags.

In [7]:
node1 = microvessel_chaste.population.vessel.VesselNode2.Create(0.0, 400.0, 0.0, reference_length)
node2 = microvessel_chaste.population.vessel.VesselNode2.Create(2000.0, 400.0, 0.0, reference_length)
node1.GetFlowProperties().SetIsInputNode(True)
node1.GetFlowProperties().SetPressure(Owen11Parameters.mpInletPressure.GetValue("User"))
node2.GetFlowProperties().SetIsOutputNode(True);
node2.GetFlowProperties().SetPressure(Owen11Parameters.mpOutletPressure.GetValue("User"))
node3 = microvessel_chaste.population.vessel.VesselNode2.Create(2000.0, 1600.0, 0.0, reference_length)
node4 = microvessel_chaste.population.vessel.VesselNode2.Create(0.0, 1600.0, 0.0, reference_length)
node3.GetFlowProperties().SetIsInputNode(True)
node3.GetFlowProperties().SetPressure(Owen11Parameters.mpInletPressure.GetValue("User"))
node4.GetFlowProperties().SetIsOutputNode(True)
node4.GetFlowProperties().SetPressure(Owen11Parameters.mpOutletPressure.GetValue("User"))
vessel1 = microvessel_chaste.population.vessel.Vessel2.Create(node1, node2)
vessel2 = microvessel_chaste.population.vessel.Vessel2.Create(node3, node4)
network = microvessel_chaste.population.vessel.VesselNetwork2.Create()
network.AddVessel(vessel1)
network.AddVessel(vessel2)

Again, we can visualize and write the network

In [8]:
scene.SetVesselNetwork(network)
scene.GetVesselNetworkActorGenerator().SetEdgeSize(20.0)
network.Write(file_handler.GetOutputDirectoryFullPath() + "initial_network.vtp")
nb_manager.vtk_show(scene, height=600, width = 1000)
Out[8]:

Next, set up the cell populations. We will setup up a population similar to that used in the Owen et al., 2011 paper. That is, a grid filled with normal cells and a tumour spheroid in the middle. We can use a generator for this purpose. The generator simply sets up the population using conventional Cell Based Chaste methods.

In [9]:
cell_population_genenerator = microvessel_chaste.population.cell.Owen11CellPopulationGenerator2()
cell_population_genenerator.SetRegularGrid(grid)
cell_population_genenerator.SetVesselNetwork(network)
tumour_radius = 300.0 * 1.e-6 * metre()
cell_population_genenerator.SetTumourRadius(tumour_radius)
cell_population = cell_population_genenerator.Update()

Again, we visualize. First turn off the grid edges so it is easier to see the cells.

In [10]:
scene.GetRegularGridActorGenerator().SetShowEdges(False)
scene.GetRegularGridActorGenerator().SetVolumeOpacity(0.0)
scene.SetCellPopulation(cell_population)
scene.GetCellPopulationActorGenerator().GetDiscreteColorTransferFunction().AddRGBPoint(1.0, 0.0, 0.0, 0.6)
scene.GetCellPopulationActorGenerator().SetPointSize(20)
scene.GetCellPopulationActorGenerator().SetColorByCellMutationState(True)
nb_manager.vtk_show(scene, height=600, width = 1000)
Out[10]:

Next set up the PDEs for oxygen and VEGF. Cells will act as discrete oxygen sinks and discrete vegf sources.

In [11]:
oxygen_pde = microvessel_chaste.pde.LinearSteadyStateDiffusionReactionPde2_2()
oxygen_pde.SetIsotropicDiffusionConstant(Owen11Parameters.mpOxygenDiffusivity.GetValue("User"))
cell_oxygen_sink = microvessel_chaste.pde.CellBasedDiscreteSource2()
cell_oxygen_sink.SetLinearInUConsumptionRatePerCell(Owen11Parameters.mpCellOxygenConsumptionRate.GetValue("User"))
oxygen_pde.AddDiscreteSource(cell_oxygen_sink)

Vessels release oxygen depending on their haematocrit levels

In [12]:
vessel_oxygen_source = microvessel_chaste.pde.VesselBasedDiscreteSource2()
#oxygen_solubility_at_stp = Secomb04Parameters.mpOxygenVolumetricSolubility.GetValue("User") * GenericParameters.mpGasConcentrationAtStp.GetValue("User")
#vessel_oxygen_concentration = oxygen_solubility_at_stp * Owen11Parameters.mpReferencePartialPressure.GetValue("User")
vessel_oxygen_concentration = 0.02768 * mole_per_metre_cubed()
vessel_oxygen_source.SetReferenceConcentration(vessel_oxygen_concentration)
vessel_oxygen_source.SetVesselPermeability(Owen11Parameters.mpVesselOxygenPermeability.GetValue("User"))
vessel_oxygen_source.SetReferenceHaematocrit(Owen11Parameters.mpInflowHaematocrit.GetValue("User"))
oxygen_pde.AddDiscreteSource(vessel_oxygen_source);

Set up a finite difference solver and pass it the pde and grid.

In [13]:
oxygen_solver = microvessel_chaste.pde.FiniteDifferenceSolver2()
oxygen_solver.SetPde(oxygen_pde)
oxygen_solver.SetLabel("oxygen")
oxygen_solver.SetGrid(grid)

The rate of VEGF release depends on the cell type and intracellular VEGF levels, so we need a more detailed type of discrete source.

In [14]:
vegf_pde = microvessel_chaste.pde.LinearSteadyStateDiffusionReactionPde2_2()
vegf_pde.SetIsotropicDiffusionConstant(Owen11Parameters.mpVegfDiffusivity.GetValue("User"))
vegf_pde.SetContinuumLinearInUTerm(-1.0 * Owen11Parameters.mpVegfDecayRate.GetValue("User"))

Set up a map for different release rates depending on cell type. Also include a threshold intracellular VEGF below which there is no release.

In [15]:
normal_and_quiescent_cell_source = microvessel_chaste.pde.CellStateDependentDiscreteSource2()
normal_and_quiescent_cell_rates = microvessel_chaste.pde.MapUnsigned_ConcentrationFlowRate()
normal_and_quiescent_cell_rate_thresholds = microvessel_chaste.pde.MapUnsigned_Concentration()
quiescent_cancer_state = microvessel_chaste.population.cell.QuiescentCancerCellMutationState()
normal_cell_state = chaste.cell_based.WildTypeCellMutationState()
normal_and_quiescent_cell_rates[normal_cell_state.GetColour()] = Owen11Parameters.mpCellVegfSecretionRate.GetValue("User")
normal_and_quiescent_cell_rate_thresholds[normal_cell_state.GetColour()] = 0.27*mole_per_metre_cubed()
normal_and_quiescent_cell_rates[quiescent_cancer_state.GetColour()] = Owen11Parameters.mpCellVegfSecretionRate.GetValue("User")
normal_and_quiescent_cell_rate_thresholds[quiescent_cancer_state.GetColour()] = 0.0*mole_per_metre_cubed()
normal_and_quiescent_cell_source.SetStateRateMap(normal_and_quiescent_cell_rates)
normal_and_quiescent_cell_source.SetLabelName("VEGF")
normal_and_quiescent_cell_source.SetStateRateThresholdMap(normal_and_quiescent_cell_rate_thresholds)
vegf_pde.AddDiscreteSource(normal_and_quiescent_cell_source)

Add a vessel related VEGF sink

In [16]:
vessel_vegf_sink = microvessel_chaste.pde.VesselBasedDiscreteSource2()
vessel_vegf_sink.SetReferenceConcentration(0.0*mole_per_metre_cubed())
vessel_vegf_sink.SetVesselPermeability(Owen11Parameters.mpVesselVegfPermeability.GetValue("User"))
vegf_pde.AddDiscreteSource(vessel_vegf_sink)

Set up a finite difference solver as before.

In [17]:
vegf_solver = microvessel_chaste.pde.FiniteDifferenceSolver2()
vegf_solver.SetPde(vegf_pde)
vegf_solver.SetLabel("VEGF_Extracellular")
vegf_solver.SetGrid(grid)

Next set up the flow problem. Assign a blood plasma viscosity to the vessels. The actual viscosity will depend on haematocrit and diameter. This solver manages growth and shrinkage of vessels in response to flow related stimuli.

In [18]:
large_vessel_radius = 25.0e-6 * metre()
network.SetSegmentRadii(large_vessel_radius)
viscosity = Owen11Parameters.mpPlasmaViscosity.GetValue("User")
network.SetSegmentViscosity(viscosity);

Set up the pre- and post flow calculators.

In [19]:
impedance_calculator = microvessel_chaste.simulation.VesselImpedanceCalculator2()
haematocrit_calculator = microvessel_chaste.simulation.ConstantHaematocritSolver2()
haematocrit_calculator.SetHaematocrit(Owen11Parameters.mpInflowHaematocrit.GetValue("User"))
wss_calculator = microvessel_chaste.simulation.WallShearStressCalculator2()
mech_stimulus_calculator = microvessel_chaste.simulation.MechanicalStimulusCalculator2()
metabolic_stim_calculator = microvessel_chaste.simulation.MetabolicStimulusCalculator2()
shrinking_stimulus_calculator = microvessel_chaste.simulation.ShrinkingStimulusCalculator2()
viscosity_calculator = microvessel_chaste.simulation.ViscosityCalculator2()

Set up and configure the structural adaptation solver.

In [20]:
structural_adaptation_solver = microvessel_chaste.simulation.StructuralAdaptationSolver2()
structural_adaptation_solver.SetTolerance(0.0001)
structural_adaptation_solver.SetMaxIterations(100)
structural_adaptation_solver.SetTimeIncrement(Owen11Parameters.mpVesselRadiusUpdateTimestep.GetValue("User"));
structural_adaptation_solver.AddPreFlowSolveCalculator(impedance_calculator)
structural_adaptation_solver.AddPostFlowSolveCalculator(haematocrit_calculator)
structural_adaptation_solver.AddPostFlowSolveCalculator(wss_calculator)
structural_adaptation_solver.AddPostFlowSolveCalculator(metabolic_stim_calculator)
structural_adaptation_solver.AddPostFlowSolveCalculator(mech_stimulus_calculator)
structural_adaptation_solver.AddPostFlowSolveCalculator(viscosity_calculator)

Set up a regression solver.

In [21]:
regression_solver = microvessel_chaste.simulation.WallShearStressBasedRegressionSolver2()

Set up an angiogenesis solver and add sprouting and migration rules.

In [22]:
angiogenesis_solver = microvessel_chaste.simulation.AngiogenesisSolver2()
sprouting_rule = microvessel_chaste.simulation.Owen2011SproutingRule2()
migration_rule = microvessel_chaste.simulation.Owen2011MigrationRule2()
angiogenesis_solver.SetMigrationRule(migration_rule)
angiogenesis_solver.SetSproutingRule(sprouting_rule)
sprouting_rule.SetDiscreteContinuumSolver(vegf_solver)
migration_rule.SetDiscreteContinuumSolver(vegf_solver)
angiogenesis_solver.SetVesselGrid(grid)
angiogenesis_solver.SetVesselNetwork(network)

The microvessel solver will manage all aspects of the vessel solve.

In [23]:
microvessel_solver = microvessel_chaste.simulation.MicrovesselSolver2()
microvessel_solver.SetVesselNetwork(network)
microvessel_solver.SetOutputFrequency(5)
microvessel_solver.AddDiscreteContinuumSolver(oxygen_solver)
microvessel_solver.AddDiscreteContinuumSolver(vegf_solver)
microvessel_solver.SetStructuralAdaptationSolver(structural_adaptation_solver)
microvessel_solver.SetRegressionSolver(regression_solver)
microvessel_solver.SetAngiogenesisSolver(angiogenesis_solver)

Set up plotting

In [24]:
scene.GetCellPopulationActorGenerator().SetColorByCellData(True)
scene.GetCellPopulationActorGenerator().SetDataLabel("oxygen")
scene_modifier = microvessel_chaste.visualization.JupyterMicrovesselSceneModifier2(nb_manager)
scene_modifier.SetVtkScene(scene)
scene_modifier.SetUpdateFrequency(2)
microvessel_solver.AddMicrovesselModifier(scene_modifier)

The microvessel solution modifier will link the vessel and cell solvers. We need to explicitly tell it which extracellular fields to update based on PDE solutions.

In [25]:
microvessel_modifier = microvessel_chaste.simulation.MicrovesselSimulationModifier2()
microvessel_modifier.SetMicrovesselSolver(microvessel_solver)
update_labels = microvessel_chaste.simulation.VecString()
update_labels.append("oxygen")
update_labels.append("VEGF_Extracellular")
microvessel_modifier.SetCellDataUpdateLabels(update_labels)

The full simulation is run as a typical Cell Based Chaste simulation

In [26]:
simulator = chaste.cell_based.OnLatticeSimulation2(cell_population)
simulator.AddSimulationModifier(microvessel_modifier)

Add a killer to remove apoptotic cells

In [27]:
apoptotic_cell_killer = chaste.cell_based.ApoptoticCellKiller2(cell_population)
simulator.AddCellKiller(apoptotic_cell_killer)

Add another modifier for updating cell cycle quantities.

In [28]:
owen11_tracking_modifier = microvessel_chaste.simulation.Owen2011TrackingModifier2()
simulator.AddSimulationModifier(owen11_tracking_modifier)

Set up the remainder of the simulation

In [29]:
simulator.SetOutputDirectory("Python/TestLatticeBasedAngiogenesisLiteratePaper")
simulator.SetSamplingTimestepMultiple(5)
simulator.SetDt(0.5)

This end time corresponds to roughly 10 minutes run-time on a desktop PC. Increase it or decrease as preferred. The end time used in Owen et al. 2011 is 4800 hours.

In [30]:
simulator.SetEndTime(20.0)

Do the solve. A sample solution is shown at the top of this test.

In [31]:
simulator.Solve()

Dump the parameters to file for inspection.

In [32]:
ParameterCollection.Instance().DumpToFile(file_handler.GetOutputDirectoryFullPath()+"parameter_collection.xml")
nb_manager.add_parameter_table(file_handler)
# Tear down the test 
chaste.cell_based.TearDownNotebookTest()
name value symbol added_by description
Generic_GasConcentrationAtStp 44.6429 m^-3 mol $$C_{stp}$$ Owen2011OxygenBasedCellCycleOdeSystem Gas concentration at STP
Owen11_BasalMetabolicStimulus 1.7 Hz $$k^0_m$$ MetabolicStimulusCalculator Basal metabolic stimulus
Owen11_CellMotilityCancer 8.33333e-15 m^2 s^-1 $$D_{cancer}$$ Owen11CaUpdateRule Maximum cell motility cancer
Owen11_CellMotilityEndothelial 1.66667e-14 m^2 s^-1 $$D_{endo}$$ Owen2011MigrationRule Maximum cell motility endothelial
Owen11_CellOxygenConsumptionRate 0.216667 Hz $$k_c^{cell}$$ User Cell oxygen consumption rate
Owen11_CellVegfProductionRate 3.33333e-05 Hz $$k_8$$ Owen2011OxygenBasedCellCycleOdeSystem Basal VEGF production rate in cell
Owen11_CellVegfSecretionRate 1.66667e-13 m^-3 s^-1 mol $$k_v^{cell}$$ User Cell vegf secretion rate
Owen11_ChemotacticSensitivity 0.333333 m^5 s^-1 mol^-1 $$\chi$$ Owen2011MigrationRule Chemotactic sensitivity
Owen11_CriticalWallShearStress 0.8 Pa $$\tau_{wall}$$ WallShearStressBasedRegressionSolver Critical wall shear stress for vessel pruning
Owen11_InflowHaematocrit 0.45 dimensionless $$H_{in}$$ User Inflow haematocrit
Owen11_InletPressure 3333.05 Pa $$P_{in}$$ User Vessel network inlet pressure$
Owen11_LatticeSpacing 4e-05 m $$\Delta x$$ User Lattice spacing
Owen11_MaxCellVegfProductionRate 0.000166667 Hz $$k_{8*}$$ Owen2011OxygenBasedCellCycleOdeSystem Max VEGF production rate in cell
Owen11_MaxTimeWithLowWallShearStress 240000 s $$T_{prune}$$ WallShearStressBasedRegressionSolver Maximum vessel survivial time with low wall shear stress
Owen11_MaximumRadius 5e-05 m $$R_{MAX}$$ RadiusCalculator Maximum possible radius
Owen11_MaximumSproutingRate 4.16667e-06 Hz $$P^{max}_{sprout}$$ Owen2011SproutingRule Maximum rate of sprouting
Owen11_MinCellCycleCancer 96000 s $$T_{min}^{cancer}$$ Owen2011OxygenBasedCellCycleOdeSystem Minimum cell cycle period cancer
Owen11_MinCellCycleNormal 180000 s $$T_{min}^{normal}$$ Owen2011OxygenBasedCellCycleOdeSystem Minimum cell cycle period normal
Owen11_MinimumRadius 1e-06 m $$R_{MIN}$$ RadiusCalculator Minimum possible radius
Owen11_OutletPressure 1999.83 Pa $$P_{out}$$ User Vessel network outlet pressure
Owen11_OxygenAtHalfMaxCycleRateCancer 186.651 Pa $$C_{\phi}^{cancer}$$ Owen2011OxygenBasedCellCycleOdeSystem Oxygen partial pressure at half max cell cycle rate cancer
Owen11_OxygenAtHalfMaxCycleRateNormal 399.966 Pa $$C_{\phi}^{normal}$$ Owen2011OxygenBasedCellCycleOdeSystem Oxygen partial pressure at half max cell cycle rate normal
Owen11_OxygenAtQuiescence 1186.57 Pa $$C^{enter}_{quiesc}$$ Owen2011OxygenBasedCellCycleModel Oxygen partial pressure at quiescence
Owen11_OxygenDiffusivity 2.41667e-09 m^2 s^-1 $$D_{c}$$ User Oxygen diffusivity
Owen11_OxygenLeaveQuiescence 1306.56 Pa $$C^{leave}_{quiesc}$$ Owen2011OxygenBasedCellCycleModel Oxygen partial pressure to leave quiescence
Owen11_OxygenTensionForHalfMaxP53Degradation 591.95 Pa $$C_{p53}$$ Owen2011OxygenBasedCellCycleOdeSystem Tissue oxygen tension for half-max p53 degradation
Owen11_OxygenTensionForHalfMaxVegfDegradation 591.95 Pa $$C_{VEGF}$$ Owen2011OxygenBasedCellCycleOdeSystem Tissue oxygen tension for half-max vegf degradation
Owen11_P53EffectOnVegfProduction -3.33333e-05 Hz $$k_{8**}$$ Owen2011OxygenBasedCellCycleOdeSystem Effect of P53 on VEGF production
Owen11_P53MaxDegradationRate 0.000166667 Hz $$k_{*7}$$ Owen2011OxygenBasedCellCycleOdeSystem Max p53 degradation rate
Owen11_P53ProductionRateConstant 3.33333e-05 Hz $$k_7$$ Owen2011OxygenBasedCellCycleOdeSystem Intracellular p53 production rate constant
Owen11_PlasmaViscosity 0.0012 m^-1 kg s^-1 $$\mu_{plasma}$$ ViscosityCalculator Blood plasma viscosity
Owen11_ReferenceFlowRateForMetabolicStimulus 6.66667e-13 m^3 s^-1 $$Q_{ref}$$ MetabolicStimulusCalculator Reference flow rate for metabolic stimulus
Owen11_SensitivityToIntravascularPressure 0.5 Hz $$k_p$$ MechanicalStimulusCalculator Shrinking to intravascaulr pressure
Owen11_ShrinkingTendency 1.7 Hz $$k_s$$ ShrinkingStimulusCalculator Shrinking tendency
Owen11_TimeDeathQuiescence 240000 s $$T_{death}$$ Owen2011OxygenBasedCellCycleModel Time for death due to sustained quiescence
Owen11_VegfConventrationAtHalfMaxProbSprouting 5e-10 m^-3 mol $$V_{sprout}$$ Owen2011SproutingRule VEGF concentration at half maximal vessel sprouting probability
Owen11_VegfDecayRate 0.000166667 Hz $$\delta_{v}$$ User Vegf decay rate
Owen11_VegfDiffusivity 1.66667e-11 m^2 s^-1 $$D_{v}$$ User Vegf diffusivity
Owen11_VegfEffectOnVegfProduction 0.04 dimensionless $$j_5$$ Owen2011OxygenBasedCellCycleOdeSystem Effect of VEGF on VEGF production
Owen11_VesselOxygenPermeability 0.001 m s^-1 $$\psi_{c}$$ User Vessel permeability to oxygen
Owen11_VesselRadiusUpdateTimestep 0.1 s $$\epsilon_t$$ User Vessel radius update timestep
Owen11_VesselVegfPermeability 1.66667e-09 m s^-1 $$\psi_{v}$$ User Vessel permeability to vegf
Secomb04_OxygenVolumetricSolubility 2.3252e-07 m kg^-1 s^2 $$\alpha_{eff}$$ Owen2011OxygenBasedCellCycleOdeSystem Oxygen solubility
In [ ]: