Skip to content

TestPythonOffLatticeAngiogenesisLiteratePaper

Jump to bottom
jmsgrogan edited this page Nov 16, 2016 · 2 revisions

This tutorial is automatically generated from the file test/python/tutorials/TestPythonOffLatticeAngiogenesisLiteratePaper.py. Note that the code is given in full at the bottom of the page.

An Off Lattice Angiogenesis Tutorial

This tutorial demonstrates functionality for modelling 3D off-lattice angiogenesis in a corneal micro pocket application, similar to that described in Connor et al. 2015.

It is a 3D simulation modelling VEGF diffusion and decay from an implanted pellet using finite element methods and lattice-free angiogenesis from a large limbal vessel towards the pellet.

Off Lattice Angiogenesis Image

The Test

import unittest
import numpy as np
import chaste.core
import chaste.cell_based
chaste.init()
import microvessel_chaste
import microvessel_chaste.geometry
import microvessel_chaste.mesh
import microvessel_chaste.population.vessel
import microvessel_chaste.pde
import microvessel_chaste.simulation
from microvessel_chaste.utility import * # bring in all units for convenience

class TestOffLatticeAngiogenesis(chaste.cell_based.AbstractCellBasedTestSuite):

    def test_fixed_outer_boundary(self):

Set up output file management.

        file_handler = chaste.core.OutputFileHandler("Python/TestOffLatticeAngiogenesisLiteratePaper")
        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.

        reference_length = 1.e-6 * metre()
        reference_time = 3600.0 * second()
        reference_concentration = 1.e-9*mole_per_metre_cubed()
        BaseUnits.Instance().SetReferenceLengthScale(reference_length)
        BaseUnits.Instance().SetReferenceTimeScale(reference_time)
        BaseUnits.Instance().SetReferenceConcentrationScale(reference_concentration)

Set up the domain representing the cornea. This is a thin hemispherical shell. We assume some symmetry to reduce computational expense.

        hemisphere_generator = microvessel_chaste.geometry.MappableGridGenerator()
        radius = 1400.0e-6*metre()
        thickness = 100.0e-6*metre()
        num_divisions_x = 10
        num_divisions_y = 10
        azimuth_angle = 1.0 * np.pi
        polar_angle = 0.5 * np.pi
        domain = hemisphere_generator.GenerateHemisphere(radius/reference_length,
                                                         thickness/reference_length,
                                                         num_divisions_x,
                                                         num_divisions_y,
                                                         azimuth_angle,
                                                         polar_angle)

Set up a vessel network, with divisions roughly every 'cell length'. Initially it is straight. We will map it onto the hemisphere.

        network_generator = microvessel_chaste.population.vessel.VesselNetworkGenerator3()
        vessel_length = np.pi * radius
        cell_length = 40.0e-6 * metre()
        origin = microvessel_chaste.mesh.DimensionalChastePoint3(0.0, 4000.0, 0.0)
        network  = network_generator.GenerateSingleVessel(vessel_length, origin, int(float(vessel_length/cell_length)) + 1, 0)

        network.GetNode(0).GetFlowProperties().SetIsInputNode(True);
        network.GetNode(0).GetFlowProperties().SetPressure(Owen11Parameters.mpInletPressure.GetValue("User"))
        network.GetNode(network.GetNumberOfNodes()-1).GetFlowProperties().SetIsOutputNode(True)
        network.GetNode(network.GetNumberOfNodes()-1).GetFlowProperties().SetPressure(Owen11Parameters.mpOutletPressure.GetValue("User"))
        nodes = network.GetNodes();
        for eachNode in nodes:

            node_azimuth_angle = float(azimuth_angle * eachNode.rGetLocation().GetLocation(reference_length)[0]*reference_length/vessel_length)
            node_polar_angle = float(polar_angle*eachNode.rGetLocation().GetLocation(reference_length)[1]*reference_length/vessel_length)
            dimless_radius = (float(radius/reference_length)+(-0.5*float(thickness/reference_length)))
            new_position = microvessel_chaste.mesh.DimensionalChastePoint3(dimless_radius * np.cos(node_azimuth_angle) * np.sin(node_polar_angle),
                                                                           dimless_radius * np.cos(node_polar_angle),
                                                                           dimless_radius * np.sin(node_azimuth_angle) * np.sin(node_polar_angle),
                                                                           reference_length)
            eachNode.SetLocation(new_position)

The initial domain and vessel network now look as follows:

Off Lattice Angiogenesis Image

In the experimental assay a pellet containing VEGF is implanted near the top of the cornea. We model this as a fixed concentration of VEGF in a cuboidal region. First set up the vegf sub domain.

        vegf_domain = microvessel_chaste.geometry.Part3()
        pellet_side_length = 300.0e-6 * metre()
        origin = microvessel_chaste.mesh.DimensionalChastePoint3(-150.0,900.0,0.0)
        vegf_domain.AddCuboid(pellet_side_length, pellet_side_length, 5.0*pellet_side_length, origin)

        vegf_domain.Write(file_handler.GetOutputDirectoryFullPath()+"initial_vegf_domain.vtp", microvessel_chaste.geometry.GeometryFormat.VTP)

Now make a finite element mesh on the cornea.

        mesh_generator = microvessel_chaste.mesh.DiscreteContinuumMeshGenerator3_3()
        mesh_generator.SetDomain(domain)

mesh_generator.SetMaxElementArea(100000.0*(units::pow<3>(1.e-6*unit::metres))); mesh_generator.Update()

        mesh = mesh_generator.GetMesh()

Set up the vegf pde

        vegf_pde = microvessel_chaste.pde.LinearSteadyStateDiffusionReactionPde3_3()
        vegf_pde.SetIsotropicDiffusionConstant(Owen11Parameters.mpVegfDiffusivity.GetValue("User"))
        vegf_pde.SetContinuumLinearInUTerm(-1.0*Owen11Parameters.mpVegfDecayRate.GetValue("User"))
        vegf_pde.SetMesh(mesh)
        vegf_pde.SetUseRegularGrid(False)
        vegf_pde.SetReferenceConcentration(1.e-9*mole_per_metre_cubed())

Add a boundary condition to fix the VEGF concentration in the vegf subdomain.

        vegf_boundary = microvessel_chaste.pde.DiscreteContinuumBoundaryCondition3()
        vegf_boundary.SetType(microvessel_chaste.pde.BoundaryConditionType.IN_PART)
        vegf_boundary.SetSource(microvessel_chaste.pde.BoundaryConditionSource.PRESCRIBED)
        vegf_boundary.SetValue(3.e-9*mole_per_metre_cubed())
        vegf_boundary.SetDomain(vegf_domain)

Set up the PDE solvers for the vegf problem. Note the scaling of the concentration to nM to avoid numerical precision problems.

        vegf_solver = microvessel_chaste.pde.FiniteElementSolver3()
        vegf_solver.SetPde(vegf_pde)
        vegf_solver.SetLabel("vegf")
        vegf_solver.SetMesh(mesh)
        vegf_solver.AddBoundaryCondition(vegf_boundary)

An example of the VEGF solution is shown here:

Off Lattice Angiogenesis Image

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

angiogenesis_solver = microvessel_chaste.simulation.AngiogenesisSolver3()

        sprouting_rule = microvessel_chaste.simulation.OffLatticeSproutingRule3()
        sprouting_rule.SetSproutingProbability(1.e6* per_second())
        migration_rule = microvessel_chaste.simulation.OffLatticeMigrationRule3()
        migration_rule.SetChemotacticStrength(0.1)
        migration_rule.SetAttractionStrength(0.5)

        sprout_velocity = (50.0e-6/(24.0*3600.0))*metre_per_second() #Secomb13
        migration_rule.SetSproutingVelocity(sprout_velocity)

        angiogenesis_solver.SetMigrationRule(migration_rule)
        angiogenesis_solver.SetSproutingRule(sprouting_rule)
        sprouting_rule.SetDiscreteContinuumSolver(vegf_solver)
        migration_rule.SetDiscreteContinuumSolver(vegf_solver)
        angiogenesis_solver.SetVesselNetwork(network)
        angiogenesis_solver.SetBoundingDomain(domain)

Set up the MicrovesselSolver which coordinates all solves. Note that for sequentially coupled PDE solves, the solution propagates in the order that the PDE solvers are added to the MicrovesselSolver.

        microvessel_solver = microvessel_chaste.simulation.MicrovesselSolver3()
        microvessel_solver.SetVesselNetwork(network)
        microvessel_solver.AddDiscreteContinuumSolver(vegf_solver)
        microvessel_solver.SetOutputFileHandler(file_handler)
        microvessel_solver.SetOutputFrequency(5)
        microvessel_solver.SetAngiogenesisSolver(angiogenesis_solver)
        microvessel_solver.SetUpdatePdeEachSolve(False)

Set the simulation time and run the solver. The result is shown at the top of the tutorial.

        chaste.cell_based.SimulationTime.Instance().SetEndTimeAndNumberOfTimeSteps(100.0, 10)
        microvessel_solver.Run()

if __name__ == '__main__':
    unittest.main()

Code

The full code is given below

File name TestPythonOffLatticeAngiogenesisLiteratePaper.py 

import unittest
import numpy as np
import chaste.core
import chaste.cell_based
chaste.init()
import microvessel_chaste
import microvessel_chaste.geometry
import microvessel_chaste.mesh
import microvessel_chaste.population.vessel
import microvessel_chaste.pde
import microvessel_chaste.simulation
from microvessel_chaste.utility import * # bring in all units for convenience

class TestOffLatticeAngiogenesis(chaste.cell_based.AbstractCellBasedTestSuite):

    def test_fixed_outer_boundary(self):

        file_handler = chaste.core.OutputFileHandler("Python/TestOffLatticeAngiogenesisLiteratePaper")
        chaste.core.RandomNumberGenerator.Instance().Reseed(12345)

        reference_length = 1.e-6 * metre()
        reference_time = 3600.0 * second()
        reference_concentration = 1.e-9*mole_per_metre_cubed()
        BaseUnits.Instance().SetReferenceLengthScale(reference_length)
        BaseUnits.Instance().SetReferenceTimeScale(reference_time)
        BaseUnits.Instance().SetReferenceConcentrationScale(reference_concentration)

        hemisphere_generator = microvessel_chaste.geometry.MappableGridGenerator()
        radius = 1400.0e-6*metre()
        thickness = 100.0e-6*metre()
        num_divisions_x = 10
        num_divisions_y = 10
        azimuth_angle = 1.0 * np.pi
        polar_angle = 0.5 * np.pi
        domain = hemisphere_generator.GenerateHemisphere(radius/reference_length,
                                                         thickness/reference_length,
                                                         num_divisions_x,
                                                         num_divisions_y,
                                                         azimuth_angle,
                                                         polar_angle)

        network_generator = microvessel_chaste.population.vessel.VesselNetworkGenerator3()
        vessel_length = np.pi * radius
        cell_length = 40.0e-6 * metre()
        origin = microvessel_chaste.mesh.DimensionalChastePoint3(0.0, 4000.0, 0.0)
        network  = network_generator.GenerateSingleVessel(vessel_length, origin, int(float(vessel_length/cell_length)) + 1, 0)

        network.GetNode(0).GetFlowProperties().SetIsInputNode(True);
        network.GetNode(0).GetFlowProperties().SetPressure(Owen11Parameters.mpInletPressure.GetValue("User"))
        network.GetNode(network.GetNumberOfNodes()-1).GetFlowProperties().SetIsOutputNode(True)
        network.GetNode(network.GetNumberOfNodes()-1).GetFlowProperties().SetPressure(Owen11Parameters.mpOutletPressure.GetValue("User"))
        nodes = network.GetNodes();
        for eachNode in nodes:

            node_azimuth_angle = float(azimuth_angle * eachNode.rGetLocation().GetLocation(reference_length)[0]*reference_length/vessel_length)
            node_polar_angle = float(polar_angle*eachNode.rGetLocation().GetLocation(reference_length)[1]*reference_length/vessel_length)
            dimless_radius = (float(radius/reference_length)+(-0.5*float(thickness/reference_length)))
            new_position = microvessel_chaste.mesh.DimensionalChastePoint3(dimless_radius * np.cos(node_azimuth_angle) * np.sin(node_polar_angle),
                                                                           dimless_radius * np.cos(node_polar_angle),
                                                                           dimless_radius * np.sin(node_azimuth_angle) * np.sin(node_polar_angle),
                                                                           reference_length)
            eachNode.SetLocation(new_position)

        vegf_domain = microvessel_chaste.geometry.Part3()
        pellet_side_length = 300.0e-6 * metre()
        origin = microvessel_chaste.mesh.DimensionalChastePoint3(-150.0,900.0,0.0)
        vegf_domain.AddCuboid(pellet_side_length, pellet_side_length, 5.0*pellet_side_length, origin)

        vegf_domain.Write(file_handler.GetOutputDirectoryFullPath()+"initial_vegf_domain.vtp", microvessel_chaste.geometry.GeometryFormat.VTP)

        mesh_generator = microvessel_chaste.mesh.DiscreteContinuumMeshGenerator3_3()
        mesh_generator.SetDomain(domain)
        mesh = mesh_generator.GetMesh()

        vegf_pde = microvessel_chaste.pde.LinearSteadyStateDiffusionReactionPde3_3()
        vegf_pde.SetIsotropicDiffusionConstant(Owen11Parameters.mpVegfDiffusivity.GetValue("User"))
        vegf_pde.SetContinuumLinearInUTerm(-1.0*Owen11Parameters.mpVegfDecayRate.GetValue("User"))
        vegf_pde.SetMesh(mesh)
        vegf_pde.SetUseRegularGrid(False)
        vegf_pde.SetReferenceConcentration(1.e-9*mole_per_metre_cubed())

        vegf_boundary = microvessel_chaste.pde.DiscreteContinuumBoundaryCondition3()
        vegf_boundary.SetType(microvessel_chaste.pde.BoundaryConditionType.IN_PART)
        vegf_boundary.SetSource(microvessel_chaste.pde.BoundaryConditionSource.PRESCRIBED)
        vegf_boundary.SetValue(3.e-9*mole_per_metre_cubed())
        vegf_boundary.SetDomain(vegf_domain)

        vegf_solver = microvessel_chaste.pde.FiniteElementSolver3()
        vegf_solver.SetPde(vegf_pde)
        vegf_solver.SetLabel("vegf")
        vegf_solver.SetMesh(mesh)
        vegf_solver.AddBoundaryCondition(vegf_boundary)

        sprouting_rule = microvessel_chaste.simulation.OffLatticeSproutingRule3()
        sprouting_rule.SetSproutingProbability(1.e6* per_second())
        migration_rule = microvessel_chaste.simulation.OffLatticeMigrationRule3()
        migration_rule.SetChemotacticStrength(0.1)
        migration_rule.SetAttractionStrength(0.5)

        sprout_velocity = (50.0e-6/(24.0*3600.0))*metre_per_second() #Secomb13
        migration_rule.SetSproutingVelocity(sprout_velocity)

        angiogenesis_solver.SetMigrationRule(migration_rule)
        angiogenesis_solver.SetSproutingRule(sprouting_rule)
        sprouting_rule.SetDiscreteContinuumSolver(vegf_solver)
        migration_rule.SetDiscreteContinuumSolver(vegf_solver)
        angiogenesis_solver.SetVesselNetwork(network)
        angiogenesis_solver.SetBoundingDomain(domain)

        microvessel_solver = microvessel_chaste.simulation.MicrovesselSolver3()
        microvessel_solver.SetVesselNetwork(network)
        microvessel_solver.AddDiscreteContinuumSolver(vegf_solver)
        microvessel_solver.SetOutputFileHandler(file_handler)
        microvessel_solver.SetOutputFrequency(5)
        microvessel_solver.SetAngiogenesisSolver(angiogenesis_solver)
        microvessel_solver.SetUpdatePdeEachSolve(False)

        chaste.cell_based.SimulationTime.Instance().SetEndTimeAndNumberOfTimeSteps(100.0, 10)
        microvessel_solver.Run()

if __name__ == '__main__':
    unittest.main()

Clone this wiki locally