Dynamics of Magnetic Nanoparticles in Newly Formed Microvascular Networks Surrounding Solid Tumours: A Parallel Programming Approach

In this paper we extend a previous 2D parallel implementation of a continuous-discrete model of tumour-induced angiogenesis. In particular, we examine the transport and capture of magnetic nanoparticles through a newly formed vascular network of blood vessels. We demonstrate how our models can be used to describe the dynamics of magnetic nanoparticles in a microvasculature and observe that the orientation of the blood vessels with respect to the magnetic force crucially affects particle capture rates leading to heterogeneous particle distributions. In addition, efficiency of magnetic particle capture depends on the ratio between the magnetic velocity and blood vessel aspect ratio. Such simulations allow a more detailed understanding of the use of magnetic nanoparticles as a mechanism for targeted anti-cancer drug delivery.


Introduction
In order to progress from the relatively harmless avascular phase to the potentially lethal vascular state, solid tumours must induce the growth of new blood vessels from existing ones, a process known as angiogenesis. To monitor and supply sufficient amounts of essential nutrients to the surrounding tissues, blood vessels have hypoxia-induced sensors, or receptors that assist in vessel remodelling to adjust the blood flow accordingly. A key mechanism of antiangiogenic therapy is to interfere with the process of blood capillary growth and literally starve the tumour of its blood supply. Indeed, a new class of cancer treatments that block angiogenesis have recently been approved and available to treat cancers of the colon, kidney, lung, breast, liver, brain, ovaries and thyroid [1][2][3][4][5].
Mathematical and computational models of vascular formation have generated a basic understanding of the processes of capillary assembly and morphogenesis during tumour development and growth [6,7]. However, by the time a tumour has grown to a size whereby it can be detected by clinical means, there is a strong likelihood that it has already reached the vascular growth phase and developed its own blood microcirculatory network. For this reason, a thorough understanding of the behavioural processes of angiogenesis is essential. Over the past few years in silico experiments focused on tumour growth have become more readily accepted by the biological community both as a means to direct new research and a route to integrate multiple experimental measurements in order to generate new hypotheses and testable predictions. This recent shift has been partly driven by the emergence of new theoretical approaches, such as hybrid modelling [8]. Hybrid models integrate both continuous and discrete processes of biological phenomena on various temporal and spatial scales. These models represent cells as individual discrete entities and often use continuous nutrient concentrations to model cellular behaviour due to their microenvironment. The cell centric nature of hybrid models naturally connects with cell biology and readily incorporates intra and extracellular phenomenon.
Recently, the use of parallel processing has highlighted the potential advantages gained from the numerical solution of complex mathematical models using high performance computing (HPC) [9][10][11][12]. HPC has evolved dramatically, in particular because of the accessibility to graphics processing units (GPUs) and the emergence of GPU-CPU heterogeneous architectures, which have led to a fundamental shift in parallel programming. Finite difference methods (FDM), such as those developed here, are the first port of call for solving complex biological phenomenon described by partial differential equations (PDEs). However, they require intensive computational resources which generally lead to significant and time-consuming expense. The advantages of time-stepping techniques in finite difference schemes lend themselves well to exploitation in a completely data parallel context. As such, parallel processing can be used to greatly accelerate such numerical simulations and offer an extremely valuable computational technique for tackling these types of problems.
In a previous paper the authors developed a 2D finite difference approximation to a hybrid continuous-discrete model of tumour-induced angiogenesis which was subsequently implemented on a parallel platform [10]. The first model presented in this paper is of a hybrid type in which a system of coupled nonlinear PDEs describe the continuous chemical and macromolecular dynamics and a discrete cellular automata-like model controls cell migration and interaction of neighbouring endothelial cells. Further, the model describes the formation of a vessel sprout network via endothelial cell migratory and proliferative responses to external chemical stimuli (i. e., tumour angiogenic factors) supplied by a nearby solid tumour, and also the endothelial cell interactions with the extracellular matrix (ECM). Once the network of blood vessels supplying the tumour has been formed, a second model based on the principles of computational fluid dynamics (CFD) is developed, and again, implemented in parallel. This model mimics the injection of drug-loaded magnetic nanoparticles (MNPs) into a primary blood vessel close to the newly formed microvascular network which is subsequently transported towards the tumour using an applied magnetic field. In this way, our second PDE model allows us to investigate the transport and capture of particles within the newly formed vasculature has a mechanism for targeted drug delivery.

Biological Description
Solid tumours generally undergo a period of avascular growth, after which they become dormant for a sustained period without access to a sufficient supply of essential nutrients, such as oxygen and glucose. Beyond a certain size (~2 mm) diffusion alone is insufficient for the provision of such nutrients; the surface area to volume ratio is too low and as such the developing tumour begins to starve. In response to this state of hypoxia, cancer cells send out signals to cells of nearby blood vessels by secreting a number of chemicals, known collectively as tumour angiogenic factors (TAF) [13][14][15]. Tumour angiogenesis stimulators include chemicals that belong to fibroblast growth factor (FGF) and vascular endothelial growth factor (VEGF) families. One important function of FGF is the promotion of endothelial cell proliferation and the physical organisation of endothelial cells into tube-like structures. Once secreted, TAF diffuse into the surrounding tissue and set up an initial steady state concentration gradient between the tumour and any preexisting vasculature. Endothelial cells situated in nearby parent vessels degrade their own basal lamina and begin migrating into the ECM [16,17]. The ECM is a complex mixture of macro-molecules, containing collagens, fibronectin etc., which functions as a scaffold for endothelial cells to grow on. The degradation of the basal lamina leads to damage, and potential rupture, of the parent vessel basement membrane. Such damage allows fibronectin from the blood to leak from the parent vessel and diffuse into the surrounding tissue [18][19][20]. Small capillary sprouts form from several endothelial cell clusters and begin to extend towards the tumour, directed by the motion of the leading endothelial cell at the sprout tip, until the finger-like capillaries reach a certain length. At this point, they tend towards each other, and form loops before fusing together in a process known as anastomoses [13,14]. Following anastomoses, the primary loops start to bud and sprout repeating the process and further extending the newly formed capillary bed. Figure 1 shows diagrammatically the general shape of the capillary sprouts and their finger-like structure. Further sprout extension occurs when some of the endothelial cells on the sprout-wall begin to proliferate. Cell division is largely confined to a region just behind the cluster of endothelial cells that constitute the sprout-tip. This process of sprout-tip migration and proliferation of sprout wall cells forms solid strands of endothelial cells within the ECM. As the sprouts approach the tumour, branching rapidly increases and produces a brush border effect, until the tumour is finally penetrated [17]. Once a supply of essential nutrients reaches the tumour, through this newly formed blood microcirculatory structure, it enters the phase of vascularisation as shown in Figure 2. To support continued growth, the vascular system constantly restructures itself implying that angiogenesis is an on-going process, continuing indefinitely until the tumour is removed or destroyed.

The Continuous Model
For a more detailed treatment of the biological aspects of tumour-induced angiogenesis as well as a more rigorous mathematical proof, readers are directed to [10,21] and references therein. Here we simply summarise the main mathematical model so as to focus on the main issues of the paper. We assume that the motion of an endothelial cell (at or near a capillary sprout tip) is influenced by three factors, namely: 1. Random motility, 2. Chemotaxis in response to TAF gradients in the surrounding connective tissue stroma, and 3. Hapotaxis in response to fibronectin gradients, also present in the surrounding tissue. So, if we denote the endothelial cell density by n, the TAF and fibronectin concentration by c and f, respectively the complete system of scaled coupled nonlinear PDEs describing tumour-induced angiogenesis can be written as: where n(x, y, t) is the endothelial cell density, D, , and are the diffusion, chemotactic, and haptotactic parameters, respectively, with c(x, y, t) and f(x, y, t) the TAF and fibronectin concentration in the 2D spatial domain (x, y) [0,1] 0,1 . All parameter values have been estimated, as far as possible, from available experimental data [21]. The system is assumed to hold on the 2D bounded spatial domain Ω (i. e. a region of tissue) with appropriate initial conditions; c(x, y, 0), f(x, y, 0) and n(x, y, 0) (see [21] for further details). Tumour cells are confined within the 2D bounded spatial domain Ω for which no-flux (Neumann) boundary conditions are imposed on Ω, the boundary of Ω, that is:

The Discrete Model
In order to capture the complex morphological features of the developing capillary network, such as individual capillary sprouts, branching and anastomosis, the continuous model must be developed further. Cellular automata models are particularly useful for providing a foundation upon which we can develop a more detailed and precise biological model. The spatial movement of individual agents in cellular automata models are primarily governed by nearestneighbour interactions and as such share some similarity with the discrete model we will present below. However, in general, the nearest-neighbour interactions for cellular automata models are based on phenomenological rules, whereas, in the discrete model presented here, the movement rules are based directly on a discretised form of the continuous model described above. The technique of tracing the path of an individual endothelial cell at a sprout tip was first proposed by Anderson et al. [22]. The method involves using standard FDM to discretise the continuous model described in (1) -(3) with the relevant boundary conditions. The resulting coefficients of the finite difference five-point stencil are used to generate the probabilities of movement of an individual endothelial cell in response to its local microenvironment. Stencil computations are those in which each node in a multi-dimensional grid is updated with a weighted average of neighbouring node values. These neighbours comprise the stencil, over which a large number of iterations across the array width generally leads to a successful numerical convergence.
We first discretise the continuous model by approximating the 2D bounded spatial domain Ω 0, 1 x 0, 1 as a grid of nodes of size h, and time t by increments of size k. By applying a forward finite difference, the fully-explicit discretised version of the continuous model (1) -(3) can be obtained. For illustrative purposes, we present the discretisation for the endothelial cell density, n: Expressions for f and c can be found in [21]. The coefficients P 0 -P 4 can be thought of as being proportional to the probabilities of endothelial cell movement i. e. the coefficient P 0 , is proportional to the probability of no movement, and the coefficients P 1 , P 2 , P 3 and P 4 , are proportional to the probabilities of moving left, right, up and down, respectively. The exact form of P 0 -P 4 can be found in [21]. We use a constant iteration size of 2,000 time steps to allow for an adequate convergence of the numerical solution. At each time step the numerical simulation involves solving Surrounding Solid Tumours: A Parallel Programming Approach the discrete model to generate the five coefficients P 0 -P 4 , which are subsequently normalised using the values of these coefficients, a set of five probability ranges are determined based on the following: . * 0 to ) * (6) where m = 1…4. A uniform random number is then generated on the interval [0, 1], and, depending on the range into which this value falls, the current individual endothelial cell will remain stationary (R o ), move left (R 1 ), right (R 2 ), move up (R 3 ), or down (R 4 ). Each endothelial cell is therefore restricted to move to one of its four orthogonal neighbouring grid nodes or remain stationary at each time step. We further assume that the motion of an individual endothelial cell located at the tip of a capillary sprout governs the motion of the whole sprout. This is not considered unreasonable since the remaining endothelial cells lining the sprout-wall are contiguous [23]. We further assume that each sprout tip has a probability, P b of generating a new sprout (branching) and that this probability is dependent on the local TAF concentration. It is also reasonable to assume that the newly formed sprouts do not branch until there is a sufficient number of endothelial cells near their tip. We will assume that the density of endothelial cells required for branching is inversely proportional to the concentration of TAF, since new sprouts become much shorter as the tumour is approached [15]. Based on these assumptions we can write down the following three cellular rules: Rule 1: New sprouts reach maturation after a length of time (ψ = 0.5) [21] before branching, Rule 2: Sufficient local space exists for a new sprout to form, and . We also assume that if a sprout tip encounters another sprout, then anastomosis can occur and a loop is formed. As a result of a tip-to-tip anastomosis, only one of the original sprouts continues to grow (purely random) and the other fuses to form the loop [24]. After the TAF has reached the parent vessel, the endothelial cells within the vessel develop into several cell clusters which eventually form sprouts [21]. For simplicity, we assume that initially five clusters develop along the x-axis at y ≈ 1, with a circular tumour located at y = 0 and the parent vessel of the endothelial cells at y = 1 as shown in Figure 3. In addition, endothelial cell doubling time was estimated at 18 hrs [25] and this is factored into our discrete model such that cell division occurs behind a sprout tip every 18 hrs. We assume that this has the effect of increasing the length of a sprout approximately one cell length every 18 hrs. Due to the inherent randomness of the discreet model, proliferation will occur asynchronously, as observed experimentally [24].

A Quick Overview of Blood Vessel Physiology
Blood flows from the heart through arteries, which branch and narrow into arterioles, and then branch further into capillaries where nutrients and wastes are exchanged. The capillaries then join and widen to become venules, which in turn widen and converge to become veins, which then return blood back to the heart through the great veins. Figure 4 shows a detailed description of the different blood vessel structures and how the network of blood vessels relate to each other.  Capillaries do not function on their own, but instead in a capillary bed, an interweaving network of capillaries supplying organs and tissues. The capillaries are the smallest of the blood vessels and are part of the microcirculation. The capillaries have a width of a single cell in diameter to aid in the fast and easy diffusion of gases, sugars and nutrients to surrounding tissues. Capillaries have no smooth muscle surrounding them and can have diameters less than that of red blood cells (erythrocytes) and may need to distort in order to pass through the capillaries. These small diameters of the capillaries provide a relatively large surface area for the exchange of gases and nutrients. Table 1 shows the different types of blood vessels, their main function and typical dimensions.

Nanoparticle Targeted Drug Therapy
Conventional chemotherapeutic agents often fail, not due to their inability to kill cancer cells, but because of their inability to distinguish cancer cells from healthy cells resulting in suboptimal efficacy combined with severe toxic side effects. The development of second generation molecularly targeted chemotherapeutic agents has emerged as one strategy to circumvent this lack of specificity. However, similar to their first-generation counterparts, many of these second-generation drugs are hydrophobic, making formulation difficult, and upon systemic administration, suffer from nonspecific bio distribution, rapid clearance and rapid degradation, in part because of their small size. For these reasons, many second-generation chemotherapeutic agents have largely failed in their quest for enhanced efficacy combined with reduced systemic toxicity [26][27][28]. In the past few decades, nanomedicine, the exploitation of the unique properties of nanoscale and nanostructured materials in medical applications, has been explored extensively as a promising strategy in the advancement of anticancer therapies with the ability to overcome many of the limitations common to chemotherapeutic agents [27][28][29]. Nanoparticles (NPs) have the potential to improve the bio distribution of chemotherapy drugs by protecting them from degradation, delivering them directly to the tumour site and preventing them from affecting healthy tissues. NPs ranging in size from 10-500 nm have been designed as drug delivery vehicles from a wide variety of materials including lipid-based amphiphiles (liposomes, hexasomes, cubosomes), metallics (iron oxide, gold), carbon nanotubes mesoporous silicates, or polymers (polymer-based micelles, drug carriers, dendrimers), as shown in Figure 5. These systems are designed such that chemotherapeutics are either physically encapsulated within or chemically conjugated to the NP.

Passive and Active Targeting
In order to specifically target nanomedicines to tumours, different approaches have been adapted, with passive and active targeting of cancer cells having been shown to be valid approaches in preclinical and clinical studies [26,28,29]. Passive targeting exploits the pathophysiological properties of the tumour vasculature which is generally highly disorganised with enlarged gap junctions between endothelial cells and compromised lymphatic drainage allowing for the extravasation of nanocarriers with sizes up to several hundred nanometres (see Figure 6(a)). Passive targeting is largely dependent on the ability of a drug nanocarrier to exhibit an increased circulation lifetime resulting in enhanced accumulation at the target site. Circulation time is dictated by the NP physicochemical properties (size, charge, biodegradability, solubility, shape, rigidity), which can be easily manipulated in the majority of the delivery systems [26,30]. The most common modification used to evade macrophage capture and increase circulation time is accomplished by making the NP surface hydrophilic through the addition of a polyethylene glycol (PEG) coating on the surface [27,28,30]. The majority of the NP-drug formulations used clinically and in development rely mainly on passive targeting.
As a means of increasing recognition of target cells by NPs, active targeting has been suggested (see Figure 6(b)). Active targeting utilises specific ligands such as peptides or antibodies that bind to molecules specifically expressed or overexpressed on target cells. Thus, active targeting does not actually improve overall accumulation at the tumour site, but rather enhances cellular uptake of the particles following their passive extravasation due to the leaky vasculature [27,28,31]. Transferrin and folate ligands are two examples of commonly used active targeting moieties in nanomedicine formulations targeting tumours [32,33]. Despite the ample evidence and extensive research effort supporting the benefits of both passively and actively targeted nanomedicines in the treatment of cancer, clinically, both strategies have met with only moderate success. This is likely due to the fact that the complexity of the tumour microenvironment (tumour heterogeneity, vascularity, location) is commonly overlooked and has a major effect on NP extravasation, accumulation, and penetration into the tumour. The tumour microenvironment is highly heterogeneous in composition with as much as half of its volume occupied by noncancerous cells and dense ECM [28]. The microenvironment creates a number of barriers that prevent these delivery systems from effectively accessing tumour cells. For example: 1. The leaky nature of the endothelium can be variable thereby restricting access to certain areas of the tumour. 2. Once NPs have exited the vessels, they usually have to pass through other cellular layers including smooth muscle cells and fibroblasts before gaining access to the tumour cells. 3. Interstitial pressure increases with increasing distance from the vessel, which can prevent NPs from penetrating deeply into the tumour. 4. Dense ECM can present an additional barrier to movement of NPs into the tumour with stiffer tumours more difficult to penetrate. The high cell density of tumour cells is difficult to penetrate, with most chemotherapy drugs only able to travel 3-5 cell diameters into the tumour and larger NPs hindered to an even greater extent.

5.
Heterogeneity in tumour cells creates challenges for active targeting as they can possess highly varied cell surface molecule expression. One strategy that has been employed which can circumvent many of the barriers encountered by NPs upon extravasation from the tumour vessels is to target NPs to the tumour vasculature. Tumour blood vessels tend to express or overexpress certain cell surface and ECM proteins that are either not present or present only at low levels in normal vessels, making them ideal as potential targets. Since the luminal surface of tumour vessels is completely accessible to circulating compounds, NPs targeting the tumour endothelium can bind to their target molecules without the need to penetrate into the tumour to deliver their contents. Recently, magnetic fields have been explored for enhancing NP delivery and efficacy in tumours [34][35][36].

Magnetic Nanoparticles
There are many different approaches to targeted drug delivery, which are classified broadly into three main categories: (i) Physical (or mechanical) approach which requires formulation of the drug using a particulate delivery device, for example a magnet which by virtue of its physical localisation will allow differential release of the drug. (ii) Biological approach which involve delivery of the drug using a carrier system like antibodies, lecithin. (iii) Chemical approach which incorporates chemical delivery systems, allow targeting of active biological molecules to specific target sites or organs, based on enzymatic activation. In magnetic drug targeting (MDT), MNPs with surfacebound drug molecules are injected into the vascular system upstream from the malignant tissue (see Figure 7), and are captured at the tumour via a localised magnetic field (usually a small rare-earth magnet e. g., neodymium magnet (NdFeB).
The NPs are manipulated precisely by applying an alternating magnetic field to transport them along the blood vessel towards the tumour. It is also possible to know exactly where they are moving with great precision by making use of the wavelength of their fluorescent emissions from the bio reactive substance they can be decorated with. The unique architecture of a tumours' blood supply makes it easy for them to absorb NPs. Instead of having a nice continuous sheet of cells as in normal blood vessels, the arrangement in tumours is very chaotic and disorganised leaving gaps. These gaps are up to 300nm, so as long as NPs are smaller than that, they will leave the blood vessel and enter the tumour. Once inside the tumour, the MNPs can be further controlled using a localised magnetic field.
Upon achieving a sufficient concentration inside the tumour, the drug molecules are released from the carriers by changing physiological conditions such as pH, osmolality, temperature, or by enzymatic activity. The released drug is taken up by the malignant cells, and the magnetic carriers are ultimately processed by the body. Since the therapeutic agents are localised to regions of diseased tissue, higher dosages can be applied which enables more effective treatment. This is in contrast to conventional therapy in which a drug is distributed in a systemic fashion throughout the body, which can have a detrimental effect on healthy tissue.
MNPs in the form of superparamagnetic iron oxide NPs (SPIONs) have received increased attention due to their characteristic small size (<10 nm) and as the name suggests, due to their superparamagnetic properties. SPIONs exist mostly as magnetite (Fe 3 O 4 ) and they can be manipulated by an external magnetic field (or magnetic field gradient) [37]. The main distinctive attribute of SPIONs is that they are superparamagnetic i. e. they generate a high magnetic moment in the presence of an external magnetic field. The remarkable response of SPION to a magnetic flux density allows for the guidance and retention of significant concentrations of the therapeutic moieties at the desired site. Furthermore, the superparamagnetic property allows these particles to convert magnetic energy to heat forming the basis of magnetic hyperthermia. Compared to other delivery methods, MNPs have a number of advantages for drug delivery because of their demonstrated responsiveness to external magnetic fields, relative safety, and versatility. Moreover, MNPs have been approved for clinical use for over a decade as magnetic resonance imaging (MRI) contrast agents and, therefore, are one of the better understood nanotechnologies in terms of patient safety. In addition, since magnetic NPs are compatible with a wide range of existing drug platforms, they can be used to effectively deliver a wide variety of therapeutic agents.

A Mathematical Magnetophoretic Model
Over the past decade, several mathematical and computational models have been developed to analyse MDT using MNPs [37][38][39][40][41][42][43][44]. Previous studies of magnetically targeted drug delivery have considered tracking individual particles under the influence of Stokes drag and a magnetic force alone. Here we also consider interactions and collisions between red blood cells within the bloodstream which cause a diffusive motion of the magnetic particles much greater than the standard Brownian diffusion. We formulate a 2D model in the same manner as [37], suitable for studying the dynamics of magnetic particles within a network of blood vessels.
Particle transport in a magnetophoretic system is governed by several factors, including: 1. Magnetic force, 2. Viscous (Stoke's) drag, 3. Inertial force, 4. Gravity, 5. Buoyancy, 6. Thermal kinetics (Brownian motion), 7. Particle-fluid interactions (perturbations to the flow field), and 8. Interparticle effects, including: I Magnetic dipole interactions, II Electric double-layer interactions, and III Van der Walls force. For most magnetophoretic applications involving submicron particles (i. e., NPs), the magnetic and viscous forces are dominant, and one can ignore all other effects. For example, the gravitational and buoyancy forces, F g and F b , respectively are given by: So, for a magnetite (Fe 3 O 4 ) particle (: 9 0.5 μm) in water ( 9 = 5,000 Kg m -3 , = = 1,000 Kgm -3 , and g = 9.81 ms -1 , we get F g = 2.57 x 10 -2 pN and F b = 0.514 x 10 -2 pN, which are more than a magnitude smaller than the applied magnetic force (~ 13 pN -see below) and subsequently can be neglected in our analysis. Similarly, the inertial force, F $ B 9 C 9 is a second order term and can therefore be neglected.

The Stoke's Drag on a Particle
The Stoke's drag, F s acting on a spherical object (particle) moving slowly through a quiescent, viscous fluid for small Reynolds numbers (Re < 1) is given by: Where µ is the dynamic (shear) viscosity of the fluid in the vessel, r p is the radius of the particle, v p and v b are the flow velocities of the particle and fluid, respectively. The dynamic viscosity of a fluid expresses its resistance to shearing flows, where adjacent layers move parallel to each other with different speeds. We assume the particles are moving through a fluid at relatively slow speeds in a laminar flow.
Laminar flow occurs at low Reynolds numbers (Re < 1), where viscous forces are dominant, and is the normal condition for blood flow throughout most of the circulatory system. It is characterised by concentric layers of blood moving in parallel down the length of a blood vessel. The highest velocity (v max ) is found in the centre of the vessel and the lowest velocity (v=0) is found along the vessel wall. The flow profile is parabolic once laminar flow is fully developed and occurs in long, straight blood vessels, under steady flow conditions (see Figure 8).
One practical implication of parabolic, laminar flow is that when flow velocity is measured using a Doppler flowmeter, the velocity represents the average velocity of a cross section of the vessel, not the maximal velocity found in the centre of the flow stream.

The Advection-Diffusion Model
The motion of magnetic particles in the blood stream is modelled as an advection-diffusion process for the particle concentration c (x, t). We consider a blood vessel in the target region, inside which magnetic microspheres are being convected (or advected). Typically, the radius of the particle, r p , is much smaller than the diameter of the blood vessel, d. This in turn is typically much smaller than the length of the magnet, l. The distance of the magnet away from the target region is usually comparable with the length of the magnet, l (see Figure 9). In addition to these forces the particle moves in response to random thermal excitations i. e. Brownian motion. This can be taken into account by considering the Brownian diffusion coefficient, D p of a particle through in a fluid (with low Reynolds number) by applying the Stokes-Einstein equation given by: where P Q is the Boltzmann's constant, T is the absolute temperature of the fluid, is the dynamic fluid viscosity and r p is the particle radius. Another diffusive mechanism that can influences the motion of particles in a blood vessel is share-induced diffusion. Blood is a highly concentrated fluid with red blood cells suspended in plasma where sheared cell-cell collisions give rise to random motions with a diffusive character. This in turn drives a diffusive motion of the plasma, causing plasma borne particles and solutes to experience shear-induced diffusion. The share-induced diffusion coefficient, D s is given by: where R D is a dimensionless coefficient dependent on red blood cell concentrations, r b is the radius of red blood cells and S is the local value of the fluid shear rate. Experimental estimates of R D for red blood cells at physiological hematocrits show a high degree of scatter, but a value of R D ≈5 x 10 -2 is representative [37]. The overall diffusion coefficient, D is given by the sum of the Brownian and shearinduced diffusivities, such that: Giving a diffusive flux vector, J D of the form: where c is the particle concentration. (14) describes Fick's first law of diffusion which postulates that the flux goes from regions of high concentration to regions of low concentration, with a magnitude that is proportional to the concentration gradient. Combining this with the advective flux vector given by: where v p is the particle velocity in the blood vessel. By applying the conservation of mass, we obtain the following advection-diffusion equation: We impose boundary conditions on c by relating the flux of particles on the boundary of the blood vessels to the evolution of the surface density on the vessel wall, X Y, Z , such that: where n is the outward unit normal at the boundary. Note that the velocity of blood flow, v b = 0 at the boundary which implies that, v 9 = v 1 .

The Magnetic Force on a Particle
The force, F 1 on a particle in a magnetic field, B is given by: Where m is the magnetic moment of the particle. Particles containing cores of magnetic material > 30nm in diameter generally have a permanent magnetic moment. The torque, T 1 = m × B causes such particles to rapidly align with the magnetic field so that the force, F 1 9 on a permanently magnetised particle is given by: However, magnetite particles of diameter < 30nm are generally superparamagnetic. In this case, m will depend on the local magnetic flux density B and is related to m through the use of a Langevin function, such that: where c • is a nonlinear saturating function given by: Therefore the force on the superparamagnetic particle, F 1 D is given by: For sufficiently weak magnetic fields, F 1 D can be approximated to the following: The magnetic field varies over the length determined by the magnet (i. e. L m = 0.01 -0.1 m).
Typical diameters of blood vessels in which the targeting takes place are much smaller than L m and thus magnetic forces across a vessel diameter can be considered approximately constant. If we consider a magnetite particle with density, 1 , magnetisation per unit mass, M, a magnetic field gradient, B G , and a magnetic volume fraction, we can calculate a typical value for the force exerted on a particle from the following: So, for a magnetite (Fe 3 O 4 ) particle (: 9 = 0.5 µm), with a density, 1 = 5.1 x 10 , kg m , , a typical magnetisation per unit mass of approximately 50 A m 2 kg -1 , the particles contain 10% magnetite by volume, and a magnetic field gradient of B G = 10 T m -1 , we get a typical force on the magnet of, F 1 = 13.4 pN.

Magnetic Particle Model
We start the network analysis by considering a 2-D model of a single small blood vessel, as shown in Figure 10. The vessel has width 2d and length L and particles enter the vessel at x = 0. The blood velocity profile is assumed to be Poiseuille flow with mean velocity U b , given by: With corresponding shear rate of the form: S ,t u |q| r i The true velocity profiles may be slightly blunted by the presence of blood cells, but this is not expected to be a significant effect. We nondimensionalise our model by scaling x with the vessel length L, y with the vessel half-width d, c with the inlet concentration c in and t with the average time taken for a particle to pass through the vessel, such that: We now drop the primes from the notation, and work only with the nondimensional variables and so the dimensionless version of (16) is given by: With shear Péclet number given by: The Péclet numbers describe the ratio of the rate of advection of a physical quantity by the flow to the rate of diffusion of the same quantity driven by an appropriate gradient. • w • ⁄ is the aspect ratio of the blood vessel, and X is the dimensionless magnetic velocity given by: We can show that )‰ ‹ ' )‰ Q for a range of particle diameters which implies that shear-induced diffusion dominates Brownian diffusion. Indeed, shear-induced diffusion plays an important role in the transport and capture of particles where the diffusive flux is strong enough to balance the advective flux due to the magnetic field. This can be written as a dimensionless parameter S, which expresses the ratio of the advective and diffusive flux across the blood vessel, given by: (31) in particular, if ' ' 1 the particle flux due to the magnetic force is negligible in comparison to the flux due to shearinduced diffusion.

The Network Flow Model
We now consider a simplified model magnetic particle capture in the newly formed network of blood vessels generated from the hybrid continuous-discrete model developed earlier in the zero diffusion limit, ' " ∞ . We ignore any transient effects and assume steady-state conditions. The relationship between the pressure drop • $ • % " between junctions' nodes i and j and the flux Qij of blood in the vessel connecting i and j is calculated assuming a Poiseuille flow profile in each vessel, such that: Where 2w $% and c $% are the width and length of the blood vessel joining nodes i and j, respectively. A schematic of the network flow model for four nodes is shown in Figure 12. At low Reynolds numbers (Re < 1), such as those found in the microcirculation, the flow deviates from the Poiseuille flow only in a small region about a blood vessel junction (shown in Figure 12). We next couple the network flow model with the magnetic particle model (27) in the zero diffusion limit, ' " ∞. In this limit the governing PDE changes type, from parabolic to hyperbolic, and so no boundary conditions are given except where the particle velocity is directed into the vessel, through the vessel wall i. e. c = 0 where v`• n = 0. The particle concentration at any point in the vessel is either 0, or the initial inlet value, c in since we are considering the zero diffusion limit and div v p = 0.
The particle flux per unit length onto the blood vessel wall, xy v`• n can be determined by integrating along the vessel wall, however, the total capture rate, ˜$ % ™ cannot exceed the inlet flux, ˜$ % xy , such that: Where n ij is the outward normal to the wall (on which capture occurs) of the vessel connecting nodes i and j. The inlet and outlet particle fluxes, ˜$ % xy and ˜$ % oe•ž , satisfy the following conservation relation: The division of particle flux at a simple vessel junction (or bifurcation) is shown in Figure 12. We can ignore the small amount of particle deposition which occurs inside the junction region and the effect of v m in this small region since, at low Reynolds numbers, the junction region is very short compared to the length of the network vessels, L ij . Consider the dividing streamline, which separates the fluid which passes into the two daughter vessels. Its position at the entrance to the junction region may be determined from the blood fluxes -*,( (flowing between nodes 0 and 1) and the flux -(, (flowing between nodes 1 and 2) and the flux -(,, (flowing between nodes 1 and 3). Using the Poiseuille flow profile in a channel of width 2d 0,1 , carrying a flux Q 0,1 , the position of the dividing streamline is calculated by choosing the flux of blood above it equal to Q 1,2 . The region of fluid containing particles at the outlet of the vessel connecting nodes 0 and 1 is computed from ˜* ,( oe•ž . Particles entering the junction region above the dividing streamline are assumed to enter the vessel connecting nodes 1 and 2, and particle below it are assumed to enter the vessel connecting nodes 1 and 3.
In this way, we determine the particle distribution without calculating the exact blood velocity field at each junction. Typical measurements for the physical parameters used in the 2D network model are shown in Table 2.

Hardware
For both the hybrid discrete-continuous model and 2D network model describing the magnetophoretic system, CUDA C algorithms were developed in Microsoft ® Visual Studio 2012 using CUDA version 7.0 and tested on an Nvidia GeForce ® GTX TM 780 GPU based on the Kepler GK110 architecture with Compute Capability 3.5. The Compute Capability describes the features of the hardware and reflects the set of instructions supported by the device as well as other specifications, such as the maximum number of threads per block and the number of registers per multiprocessor. Moreover, hardware design, number of cores, cache size, and supported arithmetic instructions are different for different versions of Compute Capability. Higher compute capability versions are supersets of lower (i. e., earlier) versions, so they are backward compatible. The operating system for both configurations was Windows 8.1.

The CUDA Programming Model
A graphics processing unit (GPU) trades off fast single thread performance and clock speed for high throughput and streaming capabilities. The GPU consists of an array of highly threaded multiprocessors (MP); each having their own individual streaming processors (SP) that share control logic and instruction cache. This keeps the available area of each SM relatively small, and therefore more SMs can be packed per die, as compared to the limited number of available CPU cores. Over the past several years, the new compute unified device architecture (CUDA) developed by Nvidia has completely revolutionised numerical computation on the GPU [45]. CUDA-driven applications run the sequential part of their workload on the CPU (the host), which is optimised for singlethreaded performance, while accelerating parallel processing on the GPU (the device). A GPU is connected to the host through a high speed bus slot, typically a PCI-express (PCI-e). The GPU has its own device memory, and is transferred between the GPU and host memories using programmed direct memory access (DMA), which operates concurrently between both the host and GPU. The device memory supports a very high memory bandwidth through the use of a wide data path e. g., a GPU with a 384-bit memory interface width allows twelve consecutive 32-bit words to be accessed from memory in a single cycle. While GPUs are frequently associated with graphics, they are also powerful arithmetic engines capable of processing thousands of lightweight threads in parallel. This capability makes an extremely viable candidate for performing highly intensive computations that exhibit high levels of data parallelism, such as the models developed here. Nowadays, modern GPUs can support up to 2,304 active threads concurrently per multiprocessor. So, for a GPU with 12 multiprocessors, this leads to more than 27,000 concurrently active threads. Threads on a CPU are generally heavyweight entities. The operating system must swap threads on and off CPU execution channels to provide multithreading capability. Surrounding Solid Tumours: A Parallel Programming Approach Context switches (i. e., when two threads are swapped) are subsequently slow and expensive. On GPUs, threads are extremely lightweight. In a typical system, thousands of threads are queued up for work in sets of 32 threads each (i. e., warps). If the GPU must wait on one warp of threads, it simply begins executing work on another. Since separate registers are allocated to all active threads, no swapping of registers or other state need occur when switching among GPU threads. Resources stay allocated to each thread until it completes its execution. In short, CPU cores are designed to minimise latency for one or two threads at a time, whereas GPUs are designed to handle a large number of concurrent, lightweight threads in order to maximise throughput. The host system and the device each have their own distinct attached physical memories. As the host and device memories are separated by the PCIe bus, data in the host memory must be communicated across the bus to the device memory. Such continual data transfer usually results in memory bottlenecks which can lead to serious performance issue when developing GPUaccelerated applications. However, such performance bottlenecks can generally minimised by making intelligent use of memory management and other optimisation techniques [11]. In essence, the CUDA programming model provides an application program interface (API) that exposes the underlying GPU architecture. In CUDA C, an instruction sequence is written into a specific function known as a kernel that can be executed on a device N times in parallel by N different CUDA threads, asynchronously. Unlike a C function call, all CUDA kernel launches are asynchronous so that control returns to the CPU immediately after the CUDA kernel is invoked [45,46].

Results and Discussion
Blood flow and particle capture results obtained from the network model are shown in Figure 13 for a typical five node simulated microvascular network based on the hybrid continuous-discrete model. This example shows a branching tree network; all the blood vessels have identical width and length. A uniform magnetic force, F` F`y s is applied to the magnetic particles. Blood and particles enter the network through a single vessel at (x, y) = (0, 0) which relates to one of the initial five clusters in the hybrid continuous-discrete model that develop along the x-axis at y ≈ 1, with a circular tumour located at y = 0 and the parent vessel of the endothelial cells at y = 1 as shown in Figure 3. The blood flows through the branches of the vascular network towards the terminal nodes, where a uniform fluid pressure, p = 0 is assumed. The division of the blood flow between the branches of the network is shown in Figure 13(a). The flow is distributed evenly over the blood vessels at each level of the network. Figure 13(b) shows the number of magnetic particles captured in each blood vessel, as a percentage of the total number of particles entering the network. Overall, about 86% of the particles are trapped in this example, and 14% are washed out of the network through the terminal vessels. Despite the even distribution of blood flow in the network the particle distribution is quite heterogeneous, with no particles being trapped in some vessels. The capture pattern depends crucially on the alignment of the blood vessels with respect to the magnetic force F m . In addition, the efficiency of magnetic particle capture also depends on the ratio between the magnetic velocity and aspect ratio of the blood vessel, X • ⁄ . Other vessel junction topologies, such as converging flow at a junction, or trifurcations can lead to more complex distributions of the magnetic particles as shown in Figure 14. These situations may be described using generalisations of the above methods and will be developed in future work. Figure 14. Additional rules for determining the particle fluxes qij may be extended to deal with trifurcations and bends where magnetic particle capture can occur on both blood vessel walls.

Conclusions
We have developed two mathematical models; one describes the development of a microvascular network of blood vessels surrounding a solid tumour. The other, a magnetophoretic system describing magnetic particle dynamics within these newly formed blood vessels. We have demonstrated that these models, when combined, can be used to describe the transport and capture of magnetic particles in a network of blood vessels and observe that the orientation of the vessels with respect to the magnetic force crucially affects particle capture rates leading to heterogeneous particle distributions. We have also shown that the efficiency of magnetic particle capture depends on the ratio between the magnetic velocity and aspect ratio of the blood vessel, X • ⁄ .
Arguably, theoretical models of angiogenesis are as diverse as their experimental counterparts, and the type and value of information they provide are equally as varied. However, the model building process itself adds value to discovery by highlighting gaps in our understanding, suggesting new experiments, providing an additional framework for hypothesis testing, and by generating new hypotheses that have quantitative basis and may not be immediately apparent from experimental results alone. With the ability to track the individual behaviours of thousands of cells, perform in silico knockout experiments that are technically infeasible experimentally, and perform high-throughput / low-cost sensitivity analyses to identify key parameters in complex systems, mathematical and computational modelling serve as additional quantitative assays that complement the available experimental data. Moreover, multi-scale models that are emerging with capabilities to integrate biological processes that span spatial and temporal scales across orders of magnitude, can be constructed in such a way that complement experimental studies and address some of these important, unanswered questions in the field of cancer. Once we have a model that is validated it is possible to more efficiently predict what should be happening in a particular experiment. With such a predictive model it is much easier and faster to perform in silico experiments to test hypotheses and predictions than running time consuming and costly laboratory experiments. More recently, the advantages of supercomputing and parallel processing techniques has highlighted the speedup, amongst other benefits, from the numerical solution of complex mathematical models of tumour dynamics such as those presented here.
The use of nanomedicines in localised drug delivery has received a lot of attention over the past couple of decades and resulted in several clinically approved formulations. These systems have been shown to have a number of advantages over conventional chemotherapeutics; however, they have not yet reached their full potential as anticancer agents. This is likely due to the fact that until more recently, features of the microvasculature and tumour microenvironment that can create barriers to effective NP delivery have been largely overlooked. With improved understanding of how the tumour microenvironment affects NP delivery and distribution within tumours, strategies can be developed to better address and overcome the shortcomings of current delivery systems. Thus, future anticancer therapies using nanomedicine can be envisioned to specifically eradicate all cancer cells within the tumour while leaving normal tissue in the body virtually untouched.
Our numerical models allow us to vary a range of parameters, including chemotactic, haptotactic and random motility coefficients to change the concentration profile surrounding the tumour and microvasculature. This will result in newly formed networks of blood vessels leading to penetration and subsequent blood flow to the tumour. More importantly, this allows us to generate multiple dynamic network structures and then model the capture and transport of magnetic nanoparticles to the tumour using our numerical simulations. In order to make full use of this virtual vascularnanoparticle system it is envisaged that more powerful algorithms will be required to cope with the vast amount of parallel simulations and subsequent 3D visualisations. As a result, we are currently investigating extending a 4D matrix operator; currently known as the Darbyshire-operator, into ndimensions to allow us to take advantage of several of its useful properties to further aid in the numerical computations. This will hopefully enable us to develop more innovative parallel algorithms and realise a goal of building a fully integrated, dynamic and real-time virtual nanolab that will greatly support current cancer research.