The supply of oxygen in sufficient quantity is vital for the correct functioning of all organs in the human body, in particular for skeletal muscle during exercise. Disease is often associated with both an inhibition of the microvascular supply capability and is thought to relate to changes in the structure of blood vessel networks. Different methods exist to investigate the influence of the microvascular structure on tissue oxygenation, varying over a range of application areas, i.e. biological in vivo and in vitro experiments, imaging and mathematical modelling. Ideally, all of these methods should be combined within the same framework in order to fully understand the processes involved. This review discusses the mathematical models of skeletal muscle oxygenation currently available that are based upon images taken of the muscle microvasculature in vivo and ex vivo. Imaging systems suitable for capturing the blood vessel networks are discussed and respective contrasting methods presented. The review further informs the association between anatomical characteristics in health and disease. With this review we give the reader a tool to understand and establish the workflow of developing an image-based model of skeletal muscle oxygenation. Finally, we give an outlook for improvements needed for measurements and imaging techniques to adequately investigate the microvascular capability for oxygen exchange.
Over the last few decades, the number of non-communicable diseases, including diabetes, heart disease, renal disease, hypertension and stroke, has increased dramatically, accounting for 60% of deaths worldwide . The development of cardio-metabolic disease is associated with changes in both macrovascular and microvascular networks in most organ systems of the body, including liver, heart and skeletal muscle. The microvasculature is characterized by vessels of diameter smaller than approximately 150 µm, whose primary function is the delivery of nutrients and oxygen to the tissue. However, if the vascular structure changes, for instance through vessel rarefication, peripheral vascular resistance may increase and/or tissue oxygen delivery may be impaired . To understand the dependence between physiological factors and disease a large body of research has been conducted using a wide range of vascular beds and disease models. To assess the importance of anatomical changes, the influence of a number of structural parameters, such as capillary density, vessel surface area and tortuosity, has been studied. Tumour networks in particular have been the subject of such investigations [3,4]. Lang et al.  found that tumour networks and healthy/normal networks differ significantly in capillary diameter, i.e. 8.0 ± 1.1 and 3.9 ± 1.1 µm, respectively. The effect of diabetes on the coronary microvasculature was studied by Jenkins et al. . In a rat model they found that compared with the microvasculature of healthy animals, diabetes resulted in microvasculature with smaller internal vessel diameters in first to third branching order from the root arteries: 194 ± 15 µm versus 267 ± 23 µm in first order, 110 ± 6 µm versus 144 ± 9 µm in second order and 73 ± 4 µm versus 110 ± 17 µm in third order.
For all research where structural hallmarks of the vasculature are assessed to study structure–function links (e.g. through image-based modelling) it is important to capture the microvascular structure within the tissue in a state as close as possible to the in vivo situation. A number of different imaging techniques can be used to capture the vascular network morphology. These differ in many aspects, such as spatial resolution and tissue penetration depth that can be achieved. Since the length scales of the relevant anatomical features of microvascular systems range from a few micrometres in two dimensions (e.g. vessel cross section) to several hundred micrometres in the third dimension (vessel length), high spatial resolution of the employed imaging technique is required as well as the ability to cover a relatively large field of view. An additional challenge in imaging the (micro)vasculature is that image contrast for biological soft tissues is usually very low and/or homogeneous in extended regions of interest for various imaging techniques. Thus, contrast agents are often applied to distinguish between features of interest . While capturing and analysing the vascular morphology ex vivo may give good indication of changes in microvascular anatomy, it provides little information of the resulting changes in vascular function, in particular capillary oxygen delivery capability, which is crucial for maintaining the tissue in a healthy state. However, studying these structure–function links experimentally in vivo introduces another layer of complexity. Depending on the imaging technique employed different limitations may arise, such as the X-ray radiation dose entailed , the accessibility of imaging facilities , a soft tissue image contrast medium that can be applied for in vivo experiments , the accessible field of view (and hence the required tissue preparation) [3,9] or limitations in spatial resolution due to animal movements. Owing to these complexities mathematical modelling is often used to inform how observed structural changes might lead to functional differences and it has frequently been employed to study blood flow and oxygen delivery to tissue. The mathematical models describing the biological processes are usually well established , but have to date been employed mainly on artificial and/or restricted geometries [11–14]. Image-based and thus realistic modelling on actual morphological data from experiments is relatively novel and has only been made possible by the recent advances in high-resolution imaging techniques and the broad availability of high-performance computing.
Image-based mathematical modelling of biological processes poses a number of challenges because various requirements have to be satisfied—at least partially—at the same time, as summarized in figure 1. In order to analyse the structure of a microvascular system, an imaging system should offer spatial resolutions and fields of view that allow a meaningful and/or relevant part of the tissue to be captured at sufficient detail. Specimen size, spatial resolution and field of view are the three factors that are important in imaging. However, for most imaging techniques, there exists an intrinsic trade-off between those factors, which cannot always be optimized independently. Another challenge is the typically low image contrast between different soft tissues (e.g. vasculature versus muscle fibres) so that identification of the different tissues and individual segmentation thereof can become difficult. Therefore, it is often necessary to draw on contrast agents, which are specific for certain tissue types in the ideal case. Finally, tissue contrasting and imaging result in certain costs, both time-wise and financially. An image-based model can also provide guidelines to determine which features need to be imaged and what spatial resolution is necessary to describe the biological process, e.g. oxygen exchange, with sufficient accuracy. Computational modelling is therefore important for determining the right balance between imaging at sufficiently high spatial resolution and investing time and financial resources inadequately, e.g. by imaging at unnecessarily high spatial resolutions or imaging at spatial resolutions that are too low and cannot provide meaningful results. The biological processes need to be modelled correctly, but too much unnecessary detail may lead to long computing times to find a numerical solution or even make the mathematical system too complex to be solved at all. Most importantly, the results provided by the mathematical model need to be validated experimentally, which—in some cases—can be used to provide (refined) parameters for mathematical modelling. Apart from the geometric information that is provided by the imaging of the system, based on which the system will be modelled , image data can supply information on material properties that feed directly into the mathematical model. The link between greyscales in computed tomography data and bone mineralization is well established , which for example can be used for modelling of structural mechanics.
When considering the problem of microvascular function different scales can be investigated, e.g. cellular and whole-organ scale . Here, we focus on the whole-organ scale by reviewing mathematical models of oxygen delivery in skeletal muscle based on ex vivo and in vivo imaging of the microvasculature to inform the association between anatomical characteristics and disease or in other words, to identify and characterize structure–function links that are important for microvascular systems. More specifically, this review discusses the different steps within the workflow for developing an image-based model of skeletal muscle oxygenation. As a result, the equations presented in this work are considered on the whole-organ scale. It may be beneficial to use homogenization approaches in this case , but these may not always hold, depending on the distribution of blood vessels in the muscle.
Skeletal muscle, which makes up almost half of a (healthy) individual's body mass , is one of the key tissues investigated in the origin and outcomes of cardio-metabolic disease. It is one of the most important tissues in terms of oxygen consumption, in particular during exercise. With that, skeletal muscle is an ideal tissue for a proof-of-concept study for image-based modelling of tissue oxygenation. The workflow presented in the current work can be directly translated to other organs, although the choice of imaging techniques and image processing methods may vary, as well as the general approach and specific parameters for computational modelling.
An image-based modelling workflow has been adopted, for example, by Cooper et al.  who created an image-based model of fluid flow in the lymph node. Measuring flow within the nodes is difficult, not least because of the small size of the nodes and their fine structure. However, through image-based modelling (on light sheet fluorescence images) Cooper et al.  were able to predict the flow pathways through the lymph node (direct from afferent to efferent lymphatics) and—for the first time—to quantify the dependency between fluid absorption/filtering through the blood vessel wall and the efferent lymphatic lumen fluid pressure. See figure 2 for a generic workflow for image-based modelling. Following the image acquisition (here by micro-computed tomography or µCT), the features of interest need to be segmented, either by hand or (semi-)automatically. For mathematical modelling of oxygen delivery in skeletal muscle, skeletal muscle tissue and blood vessels need to be segmented. The different features can now be characterized through quantitative morphometry and discretized for numerical modelling. The quantitative measures characterizing the morphology and the outcomes of the mathematical model can be correlated and an optimal structural measure for the description of the muscle oxygenation identified. The resulting numbers may then have to be validated against those that can be obtained from gold-standard imaging techniques, if non-standard imaging approaches have been used, and against biological experiments.
For a general understanding of image-based mathematical models of oxygen delivery in skeletal muscle, we briefly explain the biological background of oxygen delivery to soft tissue in §2. In §3 imaging techniques used to capture the vasculature are reviewed, along with related methods for tissue contrasting. In §4, approaches for mathematical modelling for oxygen supply in soft tissues are discussed in detail, with particular focus on image-based models. We further present validation of such mathematical models in §4.
2. Oxygen delivery to soft tissue
In this section, the biological background of oxygen delivery to soft tissues is briefly explained, which is necessary to implement a meaningful mathematical model for tissue oxygenation. A schematic of this process is provided in figure 3. This is followed by an overview of the current gold standards for the assessment of microvascular function in skeletal muscle. In this context, the importance of capturing and analysing the red blood cell distribution within the vascular system is highlighted.
2.1. Oxygen delivery to soft tissues
Blood is a suspension of red blood cells (RBCs) in an aqueous solution called plasma . The RBCs are present in a high concentration (approximately 5 × 106 mm−3), which represents about 45% of the plasma volume in a healthy adult. The ratio of RBC to the whole-blood volume is defined as the haematocrit. Oxygen is bound to the RBCs by the protein haemoglobin and one haemoglobin molecule can bind up to four oxygen molecules . Oxygen is released and diffuses down its partial pressure (PO2) gradient towards the muscle tissue. How well a solute is transferred through the capillary wall is described by the permeability of the wall to the solute, which in a capillary comprises a single layer of endothelium and its basement membrane, and the characteristics of the diffusing substrate (e.g. size, charge or lipophilicity). This process of diffusion is illustrated in figure 3. O2 exchange is very high due to the small molecular radius and because O2 is lipophilic , and it occurs across both arteriolar and capillary parts of the vascular tree [25–27].
In skeletal muscle cells the oxygen is utilized by the mitochondria (see figure 3) for adenosine triphosphate (ATP) generation to generate energy used for locomotion or exercise . Additionally, the released oxygen is bound and stored by the protein myoglobin, which is present in the sarcoplasm of skeletal muscle . Myoglobin can facilitate the oxygen supply, especially in the case of high demand, by freeing the stored oxygen into the muscle fibre for ATP generation. Each myoglobin molecule can bind one oxygen molecule . Blood flow has a significant effect on oxygen exchange between blood and muscle tissues in cases where the rate of diffusion of oxygen is significantly higher than the blood flow rate .
2.2. Assessment of microvascular oxygen exchange capability
Traditionally, the oxygen exchange capability of microvascular networks in skeletal muscle is assessed using different stereological measures on transverse and longitudinal cross sections of the muscle tissue. The most frequently used measures are given in table 1. Differences in the resulting numbers of those measures between healthy and diseased animals are investigated and assumed to be related to the disease (see examples in table 1).
The first such measure is the capillary density (CD), which is the number of capillaries per mm2 of muscle tissue, derived from a transverse cross section. This number is assumed to give a good estimate of the oxygen supply capability and hence aerobic metabolism of a muscle . However, the capillary density depends strongly on the diameter of muscle fibres. Namely, a muscle with larger fibres has a lower CD . The capillary-to-fibre ratio (C : F) is the mean number of capillaries that are adjacent to a single muscle fibre. Similar to the capillary density, this value is assumed to provide information on the metabolic activity of the muscle; it does however depend on the muscle fibre type. For mouse soleus muscle, Poole et al.  found a capillary-to-fibre ratio of 2.17 ± 0.06, Hudlicka  indicated 2.3 ± 0.03 and Dapp et al.  put forward a capillary-to-fibre ratio of 1.83 ± 0.04. This shows that the reported capillary-to-fibre ratio is relatively consistent. As the capillary-to-fibre ratio does not depend on muscle fibre area it may not immediately yield information about oxygen supply capability. Other measures are length and volume density, which in contrast with CD and C : F consider the three-dimensional (3D) structure of the blood vessels. They describe the length or volume of capillaries over the volume of corresponding muscle tissue. Kondo et al.  reported that rats with type 2 diabetes have a different capillary volume in the soleus muscle when compared with healthy animals (47% lower in diabetic rats). At the same time, no significant difference in CD was found.
Vessel tortuosity is a measure that is often considered when studying vasculature of organs and muscles, as it is assumed that higher tortuosity of vessels—through increase of surface area and reduction of oxygen diffusion distances—can be associated with an increased oxygen supply . There are a number of different methods used for assessing tortuosity, the most common being the vessel segment length divided by the Euclidian distance between the segment's endpoints . However, there is also a hypothesis that tortuosity is not influenced by oxygen demand and supply, but is related to the contractile state of the muscle and therefore the sarcomere length [37–39], which may in turn improve oxygen supply to the muscle tissue. Poole et al.  found that tortuosity in capillaries increases when the sarcomere length is reduced below 2.0–2.4 nm. In the same study, it has been shown that the length of the capillary network depends on the tortuosity to an extent of 24–38% . Moreover, Mathieu-Costello et al.  and Poole et al.  found no significant link between capillary tortuosity and body size, training, athletic ability or aerobic capacity. Based on these results one can conclude that tortuosity cannot be the sole indicator for differences between oxygen exchange capabilities of the microvasculature in skeletal muscle in health and disease.
All the above measures take into account only the blood vessels that have been visualized by a specific imaging technique. Differences in numbers of these measures between healthy and diseased subjects are then used to associate the disease with the capillary oxygen exchange capability. Poole et al.  oppose this approach to examine O2 supply and reviewed the literature to show that impairment of vessel function has to do with differences in flux, i.e. RBC flux, blood–myocyte oxygen flux, and haematocrit, rather than number and volume of perfused capillaries. RBC flux is present in most capillaries at rest and increases with exercise to meet the muscle demand . In agreement with this, Maeda et al.  reported that the number of capillaries exchanging oxygen was 93.0 ± 5.5% of all blood-perfused capillaries in rat soleus muscle. A higher difference in these numbers was found by Fraser et al.  in rat extensor digitorum longus muscle (EDL), i.e. around 12% of visible capillaries stayed unperfused during an observation period of 60 s.
Overall, it remains inconclusive whether or not structural measures of muscle capillarization can be used as explicit indicators of oxygen exchange capability. Mathematical modelling can help clarify this question. An association between metabolic syndrome and skeletal muscle microvascular haematocrit has also been found , thus suggesting that for investigating the impairment of microvascular oxygen exchange capability, not only the morphology of the blood vessels has to be taken into account, but also the distribution of RBCs within the vessels. Therefore, in order to assess microvascular function meaningfully, the imaging technique adopted should be able to resolve the microvascular network structure, and simultaneously the distribution of RBCs within the blood vessels, and if possible, record changes of the network structure and the RBC distribution over time. By using mathematical modelling it is possible to study the muscle tissue oxygenation of different capillary networks (i.e. tissue function) and to link the results to the structural measures of those networks (i.e. tissue structure). This will increase the understanding of structure–function relationships for (micro)vascular systems.
3. Imaging systems for vascular networks
An intrinsic problem of the previously presented measures as indicators for microvascular exchange capability is that they have only been derived for sections of muscles or small tissue volumes obtained by light microscopy and confocal laser scanning microscopy (CLSM, a technique based on laser excitation and fluorescent emittance by the sample, see below), respectively. It is unclear whether results of small sampling volumes can be extended to the whole muscle or how local differences in microvasculature are taken into account. Using CLSM, Cebasek et al.  reported that they observed significant differences between results when a two-dimensional (2D) or 3D imaging technique is employed to study rat soleus vascularization. Applying a 3D measure, i.e. length of capillaries over the length of the muscle fibre they surround, they found a capillarization up to 40% higher than that resulting from C : F measurements. Therefore, there is a need for high-resolution 3D imaging techniques to assess the microvasculature in health and disease, taking into account a representative portion of the muscle. These different requirements are difficult to meet simultaneously for most imaging systems. Below, we review imaging techniques for the assessment of microvasculature and we discuss their advantages and disadvantages, including relevant tissue contrasting methods.
Different imaging techniques consist of the same two fundamental components: a source and a detector. The source generates an excitation signal. This signal interacts with the sample lying in the signal path. The interaction of the signal with the sample either changes the original signal and/or generates a new signal emitted by the sample. The detector captures this changed and/or new signal. Excitation signals include visible light, electron beams, X-rays, magnetic pulses, positrons or acoustic waves . Different imaging techniques can be distinguished by their different characteristics, such as 2D or 3D readouts, excitation mechanism, capability to penetrate tissue and compatibility for imaging in vivo or ex vivo. Imaging systems for the visualization of vascular networks have been reviewed extensively, e.g. in the context of blood flow modelling by Kim et al. , angiogenesis by Kiessling et al.  and systems biology by Kherlopian et al. . Here, we have created a scatter graph that compares maximum spatial resolution and tissue penetration depth (figure 4). In table 2, an overview is provided of the capabilities of the different imaging techniques. These collectively show that while light microscopy is used as a gold standard for the visualization of the microvasculature (in 2D), the most suitable imaging techniques for capturing the structure of the (micro)vascular tree are µCT, synchrotron radiation-based CT (SR CT), and possibly light sheet fluorescence microscopy (LSFM). When attempting to capture the whole vascular network on an organ level, neither light microscopy nor CLSM is suitable due to their small tissue penetration depth (approx. 100 µm). For these techniques, the tissue usually needs to be physically sectioned, which is time-consuming and often leads to image artefacts due to sample distortion and deformation, caused mainly by the shear stresses resulting from cutting. If the microvasculature within a whole murine skeletal muscle (with typical dimensions of about 0.2 × 0.2 × 0.8 cm3) is to be imaged non-destructively, the imaging technique must provide a sufficient depth of field to penetrate the whole soft tissue sample. This limit of 0.2 cm is marked by a blue line in figure 4. Furthermore, only techniques with a spatial resolution of at least about 1–1.5 µm can be considered for the assessment of microvascular systems. This is the limit when the finest capillaries of 5 µm diameter will be represented by 3–5 pixels/voxels, which is necessary to have confidence when identifying and segmenting structures (due to partial volume effects) as well as for performing image-based modelling subsequently in a reliable manner, as recommended for finite-element modelling of trabecular bone microarchitecture  for instance. The dashed red line in figure 4 marks this limit.
3.1. Light and fluorescence microscopy
Light microscopy (LM) was the first method for biological imaging, introduced around the seventeenth century as a single lens microscope. Since then LM has evolved and reached its interim maximum spatial resolution of roughly 0.2 µm in the nineteenth century, which is the diffraction limit at half of the minimal wavelength of visible light. Recently, methods have been developed to image beyond this diffraction limit, at so-called super-resolution, i.e. in the range of tens of nanometres . Different forms of LM exist, such as polarized, phase-contrast or fluorescence microscopy. These imaging modalities make use of different physical principles to make features of interest visible within the material or tissue to be examined. Fluorescence microscopy is based on autofluorescence or fluorescence of a specific tissue stain to highlight features in different colours, whereas phase-contrast microscopy detects different local phase shifts of the light wave, induced by different local material properties of the specimen. Similarly, polarized microscopy detects the change in light polarization . LM is the gold standard for most imaging techniques. It has been used extensively in imaging skeletal muscle to define changes in capillarization due to disease and training [29, ch. 5], as reviewed by .
CLSM is fluorescence microscopy with the advantage that it has a greater depth of field due to the reduction of out-of-focus light via a pinhole. The depth of field for CLSM is limited to about 150 µm, due to a limitation of tissue penetration of light, while the maximal spatial resolution is roughly 100 nm . Some studies have reported a penetration depth of up to 1500 µm if the tissue is cleared chemically prior to imaging [51,52]. CLSM has been used extensively to image the 3D microvasculature in skeletal muscle [31,44,53–55]. Samples need to be autofluorescent or stained with fluorescent markers to be excited and to emit a response signal.
A relatively recent development of LM is the so-called light sheet fluorescence microscopy (LSFM), where a laser is used to create a light sheet, which excites a several micrometres thick plane within the sample. The fluorescence response of the sample is then detected . In contrast with CLSM it has the advantage of little photobleaching (fading of the fluorescent dye due to non-specific illumination), as the fluorescent response is activated only as each plane is illuminated selectively . LSFM can be employed to capture sample volumes of up to 1 cm3, without the need for any physical sectioning. The in-plane resolution of LSFM is about 1 µm . However, the technique requires the tissue to be cleared chemically for the laser sheet to be able to penetrate it, and it depends on the diffusion capability of fluorescent stains to permeate the tissue. Owing to these limitations and its general set-up in terms of size and distances of the components, LSFM has had little application in imaging organs of larger animals. Most studies investigated zebrafish and fruit fly , few looked at brain neuronal activity in the mouse [58,59] and Mayer et al.  used it to study high-endothelial venules in mouse lymph node. However, to the best of knowledge of the authors, LSFM has yet to be trialled in mouse skeletal muscle for visualization of vascularization. Combined with µCT, LSFM has the potential to yield more conclusive results as it can visualize features not normally observed in µCT, e.g. the lymphatic network, nerves and different cellular structures.
3.2. Micro-computed tomography
µCT is a non-destructive 3D imaging technique that makes use of variable X-ray absorptions induced by the different components inside a tissue. 2D radiographies or projections of a sample are collected for different angular position (but for at least 180°), which are then reconstructed in 3D to form a volumetric representation of the sample .
The main advantages of µCT are its non-destructive nature and the high spatial resolution that can be obtained . Typical maximal spatial resolutions of commercial µCT systems are just below 1 µm, where synchrotron-based CT set-ups normally operate down to the diffraction limit of visible light. One general drawback of µCT is that at high spatial resolutions the exposure to significant levels of X-ray radiation is a crucial factor to consider. Hence, the maximal pixel size is typically limited to around 5 µm for in vivo applications, resulting in actual spatial resolutions in 3D around 10 µm. Moreover, animal or sample movement during scanning detrimentally affects image quality of the reconstructed CT data, which becomes critical when high spatial resolutions are targeted. Commercial scanners exist for small animals such as mice, rats or rabbits, where the heart rate can be monitored for respiratory gating under anaesthetic, in order to avoid or at least minimize movement artefacts. Badea et al.  and more recently Clark & Badea  have reviewed the suitability of µCT for in vivo imaging of small animals. A lethal dose for mice (50% of the population die within 30 days) is estimated to be 5–7.6 Gy . For low-resolution µCT scans (e.g. 135 µm pixel size) this does not pose a problem as the reported dose at this level is around 0.25 Gy. By contrast, for high-resolution scans required to assess microvasculature (5 µm pixel size) the lethal dose can easily be reached . However, rodents can repair damages induced by radiation and die not before a cumulated dose higher than 30 Gy, a threshold six times higher than the cumulated 5 Gy after a 6-week in vivo imaging study . Furthermore, modern detector technology provides more sensitive scintillators and thus higher photon flux that is recorded by the detector, which results in shorter scanning times. This way, it may be possible to reduce the X-ray radiation dose if the detector pixel size can be further reduced for high-resolution imaging . Also, further work on optimizing tube current and voltage may help to decrease dose . In conclusion, µCT has reached a point where it is commonly used for in vivo imaging of small rodents.
µCT has been used in many studies on vascular networks; some examples include imaging of the ocular and renal microvasculature [69–71], imaging of scaffold–bone interaction , phenotyping of cardiovascular development in mouse embryos , imaging of remodelled bioactive glass foam scaffolds , functional imaging of rat hearts , studying cavernous haemangioma of the liver , examining neovascularization in tissue-engineered bone constructs  or computer-aided design of microvasculature [78,79]. Heinzer et al.  used a combination of both µCT and SR CT to perform a multi-scale analysis of vascular networks. Similarly, Schneider et al.  studied the murine hind limb vasculature using corrosion casts with both µCT and SR CT. Razavi et al.  imaged the pulmonary and cardiac circulation using µCT. A review of µCT imaging of vascular networks can be found in Bentley et al.  and more recently in Zagorchev et al. .
As imaging contrast of organic soft tissue is typically limited for all imaging techniques, it is conventionally considered necessary to prepare tissues with staining agents. The choice of staining agent strongly depends on the imaging technique to be used and the feature of interest to be assessed. It is necessary to consider several aspects, such as costs, toxicity, tissue shrinkage and perfusion range. The methods to enhance vascular contrast for µCT imaging can be differentiated into methods that stain vessel walls (§3.2.1) and the lumen of the vessel (§3.2.2).
3.2.1. Staining methods
Metscher  and Pauwels et al.  have reviewed a number of different staining agents for soft tissue. The most common stain used is osmium tetroxide (OsO4), which offers good imaging contrast for electron microscopy and µCT and induces low tissue shrinkage. OsO4 is expensive, very toxic and requires specific disposal, which results in high financial costs [85,86]. Another frequently used agent is Lugol's solution, which is based on iodine. Lugol's solution yields low toxicity, leads to low tissue shrinkage and is easy to prepare [73,85]. However, it is prone to leakage into the interstitial space. Phosphotungstic acid (PTA) is another agent which provides high contrast and stability and has low toxicity. PTA binds to collagen, fibrin and fibres of connective tissue [6,85]. Pauwels et al.  compared another 27 staining agents, among others mercury(II) chloride (HgCl2), phosphomolybdic acid, PTA and ammonium orthomolybdate ((NH4)2MoO4), concluding that most of these were similarly suitable for the visualization of soft tissue. However, depending on the agent, sample size or staining times need to be adjusted . These methods have also partly been reviewed for soft tissue staining . Overall staining methods are rather unsuitable for imaging the microvasculature in skeletal muscle, as they are largely non-specific, leaving only the vessel lumen unstained. Thus, even at high spatial resolution segmentation of capillaries becomes difficult as the vessel lumen is small and the vessels may easily collapse, making them difficult or even impossible to detect.
3.2.2. Perfusion methods
Two different sets of perfusion methods exist, firstly perfusion with polymers for corrosion casting and secondly perfusion with nanoparticles. The most commonly used polymer for corrosion casting for µCT is Mercox, commercially available in its non-hazardous form Mercox II (Ladd Research, Williston, Vermont, USA). It is used for classical corrosion casting, i.e. the surrounding organic tissue is macerated after perfusion and curing [88–90]. Mercox polymerizes quickly, which means that perfusion times and volumes must be relatively small . Mercox has good permeability properties within the entire vascular network and shows minimal shrinkage . However, one of its drawbacks is its low X-ray absorption. In order to enhance X-ray absorption, it can be coated with a material of high atomic number, such as osmium tetroxide . Microfil (Flow Tech Inc., USA) is a frequently used silicone rubber; a polymer that is fluid during perfusion and hardens within 30 min [71,82,92–97]. A problem often encountered with this agent is that it does not always perfuse all of the microvasculature . Contrary to that, Ghanavati et al.  reported a very uniform filling of the cerebro-vasculature using Microfil. Apart from Microfil, other casting agents are common, such as the polyurethane-based resin PU4ii (vasQtec, Zurich, Switzerland), for which good perfusion of the cerebral and the murine hind limb microvasculature was reported [81,99–101]. It is also possible to perfuse the vasculature with nanoparticles, which are small particles of a high atomic number material with dimensions of about 100 nm. This can be, for example, a suspension of gold nanoparticles [102,103], liposomal iodine nanoparticles  or bismuth sulfide (Bi2S3) nanoparticles. Nanoparticles can also be used for in vivo µCT imaging . The nanoparticles used for µCT require long half-lives, i.e. the time before the half of the particles have left the body/tissue of interest, because µCT scans generally take at least a few minutes and up to a couple of hours, depending on the set-up and experimental settings . The size of the nanoparticles used can be varied depending on the vasculature investigated. For example, as tumour vasculature is leaky, larger particles may be beneficial or by contrast, smaller particles can be used in order to determine tumour permeability . Owing to the segmentation difficulties arising from staining methods, perfusion methods should be the method of choice for absorption-based µCT of the microvasculature in skeletal muscle. However, as perfusion with polymers has proven unreliable, nanoparticles are a viable alternative.
3.3. Phase-enhanced imaging
Another method to gain image contrast using µCT is by exploiting the fact that X-rays are not only partially absorbed, but that the X-ray waves are also refracted at the interfaces between different sample features. This can be seen from the complex refractive index, which is given by where δ is the refractive index decrement, related to the phase shift. β is related to the attenuation coefficient μ via the X-ray wavelength λ:
For low-absorbing materials the ratio of the phase shift in comparison to the X-ray absorption is very high, e.g. δ/β = 100 for calcium in bone, whereas it is δ/β = 104 for water, e.g. in soft tissue, at an X-ray energy of 14 keV. The ratio of δ to β increases further with increasing photon energy . With regard to applications in a clinical setting, it becomes clear that phase-contrast X-ray CT has a distinct advantage over traditional CT as the ratio between phase shift and absorption increases as the X-ray energy increases, resulting in lower dose and better soft tissue contrast at the same time.
The simplest phase-contrast CT technique in terms of implementation is free space propagation, where the X-rays are allowed to propagate towards the detector in free space. This makes use of the Fresnel edge diffraction of different rays of the beam when they are hitting boundaries between different materials (see figure 5). Propagation-based phase contrast is especially useful for high-resolution phase-contrast imaging, as most other phase-contrast techniques cannot reach a resolution below a few micrometres . As most other phase-contrasting techniques, free-phase propagation is mainly used at synchrotron sources, as most laboratory-based X-ray CT sources cannot provide sufficient beam coherence. Coherence of the X-ray beam is important for the interference and hence phase information to be extracted correctly . Another advantage of synchrotron radiation includes the fact that the X-ray beam is monochromatic and is delivered at a high flux, both of which are not essential for propagation-based phase contrast imaging , but avoid beam hardening effects and reduce scanning time significantly, respectively. Phase-contrast µCT has successfully been used to image sub-micrometre particles in unstained lung tissue with a voxel size of 370 nm .
However, the technique has been used relatively little to image the vasculature. Lang et al.  have been able to identify tumour vasculature and only recently it was reported that phase-contrast imaging in conjunction with phase retrieval at synchrotron facilities enabled simultaneous visualization of blood vessels and nerve fibres in the spinal cord without the need for contrast agents . Using state-of-the-art laboratory-based µCT systems that are partially coherent due to their small spot size, Walton et al.  demonstrated that phase enhancement can also be exploited to some extent using commercial µCT scanners. Namely, Walton et al.  visualized sub-micrometre structures within skin and rat artery walls using µCT and they reported that subsequent histological and immunohistochemical staining was compatible with prior X-ray exposure.
Other phase-contrast techniques include grating interferometry, for which the requirement of coherence is lifted as the coherence can be created by placing a grating between the source and the sample . Grating interferometry makes use of two further gratings placed in between sample and detector to extract phase-shift information. Using this technique it was possible to obtain images of rat brain , images of rat testicle at a comparable quality of histological images  and differential phase and dark field images of a cancerous mastectomy , showing the potential of grating-based interferometry for pre-clinical applications. All of the mentioned experiments have been carried out ex vivo. Other experimental set-ups for grating-based interferometry and other phase-contrast imaging techniques have been reviewed by Bravin et al. . For ex vivo imaging with the purpose of providing image datasets for mathematical modelling on the whole-organ scale, phase-contrast SR CT can be considered as an optimal technique if access to a synchrotron facility can be obtained, due to the easy experimental set-up and simple specimen preparation requirements (fixation and potentially embedding). However, if one is limited to laboratory equipment, X-ray absorption-based µCT should be considered, as scanning times required for phase-contrast imaging can be lengthy. This is due to significant loss of photon flux, caused by extended sample-to-detector distances that are required for free-space propagation and by the absorption of the source grating.
3.4. Image processing for image-based modelling
As highlighted in figure 2, there are a number of steps involved before images can be used for modelling of processes that occur. In particular, the structures which are relevant for the biological processes that occur need to be identified, i.e. in the case of muscle tissue oxygenation the muscle tissue itself, as well as the features that share interfaces with the structure in a way that is significant for the biological process. In the case of tissue oxygenation these are most importantly the blood vessels which are supplying the oxygen. That is, the flow of oxygen across the boundary between the blood vessels and the muscle tissue needs to be described by the mathematical model (see (4) in figure 3) and therefore the boundary needs to be identified. This identification and extraction of structures in the images is called segmentation. At the end of the process of segmentation an adapted image (or image stack) is obtained that contains only the structures of interest, which can be distinguished by their different greyscales in the case of CT.
For image segmentation a range of different software packages are available depending on the specific requirements and characteristics of the images. For instance, these are Fiji/ImageJ , which is freeware, and a range of commercial image processing softwares, such as Avizo (FEI, USA), VGStudio (Volume Graphics GmbH, Germany) or ScanIP (Simpleware, Synopsis Inc., USA). Prior to segmentation it is advantageous to pre-process the images to facilitate segmentation. In particular, filters of different kinds can be applied to reduce noise, sharpen the images, highlight edges or remove artefacts that can vary between different imaging techniques. For example, it is often necessary to apply a ring removal algorithm after µCT/SR CT imaging, where ring artefacts are due to defect or wrongly calibrated detector elements. For an overview of image filters and automated feature segmentation approaches, see Nixon and Aguado . After pre-processing and depending on the quality of the images as well as on the image contrast the segmentation of features of interest can be performed automatically, semi-automatically or by hand, for instance by using greyscale thresholding or a region-growing algorithm (possible for blood vessels perfused with a polymer for corrosion casting) or by a combination of manual segmentation and automated methods [118,119]. For manual segmentation, the expenditure of time to segment blood vessels in an image stack containing a whole muscle, recorded at micrometre resolution, can amount to several days or even weeks. In contrast with this, region-growing segmentation will require only a number of hours, depending on the computing power, while segmentation based on global thresholding is a matter of seconds. The accuracy of the segmentation results is largely dependent on user expertise in the case of semi-automatic and manual segmentation.
Based on the segmented blood vessel network the structural parameters commonly used for characterization of microvascular oxygen exchange capability, e.g. capillary density or vessel tortuosity, can be computed. This task can often be performed using the same software package that has been used for segmentation. For ImageJ in particular a plugin has been developed for characterizing bone structure, called BoneJ , which can also be used to study the morphology of the microvasculature. In some cases, filters or algorithms need to be applied prior to the actual computation, to determine capillary number and tortuosity for instance. Namely, a skeletonization algorithm can be used (available for example in Fiji ), which reduces a tubular structure to its geometric centre line. The resulting points of the skeleton per image and muscle slice can then be counted, for example, to determine the capillary density through voxel counting or the tortuosity by dividing the Euclidian length and actual length of a vessel segment's skeleton, which are both outputs provided by the skeletonization plugin within Fiji . While being able to provide structural information, images of muscles and their microvasculature do not immediately reveal information of muscle functionality, in particular of the vasculature's ability to supply the muscle with sufficient oxygen throughout. Morphometric measures, such as capillary density, give a first insight into this, but it is only with mathematical modelling of oxygen perfusion that it becomes possible to relate such structural information to the actual ability of the microvascular network to meet the tissue oxygen demand.
To this end, it is necessary to discretize all extracted features after segmentation. Clearly, the voxelized representation of the data is one form of discretization; however, it is often preferred to create a smoother representation of the geometry. Discretization approaches include the finite element [119,121], finite cell , finite volume  or a finite difference method. One can differentiate also into approaches that discretize by defining image voxels as finite elements and by creating a mesh that approximates the realistic shape of the vascular structure with smaller elements . For 3D datasets the finite-element method is most commonly used. In this case a polygonal volume mesh of the data is created, e.g. a tetrahedral representation. Smoothing of the image data prior to meshing is sometimes required to create a suitable mesh with sufficiently large and no intersecting elements. For example, single voxels need to be removed by despeckling (e.g. by component labelling) or boundaries have to be smoothed (e.g. by applying dilation and erosion operations). Many of the commercial software packages introduced above include modules capable of generating surface or volume meshes, the former of which can be imported into many modelling packages and re-meshed into a volume mesh (OpenFOAM, OpenCFD Ltd, ESI Group, Paris, France; COMSOL Multiphysics, Comsol Ltd, Cambridge, UK). Freeware meshing software is also available, such as the tool MeshLab. In some cases, it is also possible that the software package used for mathematical modelling allows for meshing an imported geometry (COMSOL Multiphysics, Comsol Ltd, Cambridge, UK; Ansys, Ansys Inc., Canonsburg, PA, USA; Abaqus, Dassault Systèmes SE, Vélizy-Villacoublay, France). Meshing of specific and particularly complex biological geometries is a challenging task, which often requires large computational power and memory.
Another important aspect involved in image-based modelling of biological processes is the formatting of the image data. Images are often provided in a specific file format, depending on the image acquisition system, such as raw, TIFF or DICOM, but also file formats that are proprietary to the software and/or hardware of the imaging device manufacturer. To import these original images into specific image processing and/or meshing software, they might need reformatting. It is advisable to keep to a universal, non-compressed file format, such as raw or TIFF, as datasets of this type will usually be compatible with most software packages. In case that a large contrast resolution is not needed for segmentation, for example, when heavily and consistently stained tissue is distinguished from a homogeneous background, the dynamic range of the images can be reduced to save data storage and reduce network load and computing time for image pre-processing and segmentation. Depending on the meshing technique and software the output file format can vary, but most commercial software packages offer a range of output file formats, which are compatible with most of the (commercially) available modelling software packages. It is important that the dataset returned by the meshing software contains all domains and boundaries that are necessary for subsequent computational modelling. In addition to the formatting of the data, the availability of image and modelling data to the greater research community needs to be considered. It is important to make the data freely available to others so that they can reproduce results and build upon them. Furthermore, if research groups with a focus on experimental and imaging expertise provide imaging data for other groups experienced in mathematical modelling, we believe that additional knowledge can be generated, which is out of reach for both individual groups alone. Moreover, research councils are increasingly asking to make experimental data available as a condition of future funding. However, image data and output of mathematical models, especially in a 3D context, can exceed several gigabytes and terabytes in size, which can be an obstacle for simple upload to a storage server. The Materials Data Centre (MDC) funded by the UK Government (JISC)  and the German Helmholtz Association project ‘Large Scale Data Management and Analysis’ (LSDMA) [125–127] are addressing the problem known as ‘Big Data’ and propose solutions suiting different research communities.
4. Mathematical modelling of tissue oxygen supply
We review here models of the oxygen supply to soft tissues with a particular focus on image-based models. The equations and parameters needed for modelling of tissue oxygenation are given in appendix A. While the oxygen delivery to skeletal muscle does partly depend on blood flow, models of blood flow have been reviewed previously in  and to some extent in  in the context of oxygen transport and regulation. The main outputs of blood flow modelling, used for further modelling of tissue oxygen perfusion, are the partial pressure of oxygen in blood and its oxygen saturation. Partial oxygen pressure and oxygen saturation are directly related to the blood haematocrit. Most models of molecule transport in and across blood vessels use parameters and equations reviewed in the Handbook of Physiology . The theory of oxygen transport to tissue was also presented and reviewed by Popel . Jain [128,129] reviewed the transport of molecules (including oxygen) in the tumour interstitium and across the tumour vasculature, while  reviewed theoretical models of oxygen transport to tissue. In the review of  models of tumour oxygenation are discussed, which are similar to oxygenation of other soft tissues. This review focuses on the whole-organ scale of muscle tissue oxygenation and is thus only describing models that are relevant to this context.
4.1. Modelling of artificial networks
Many models of oxygen diffusion in soft tissue are implemented on artificial vascular networks, which are based on observations by imaging or theoretical approaches. The first model of oxygen transport from the capillaries to muscle tissue was developed nearly hundred years ago by Krogh in 1919. Krogh determined diffusion coefficients of oxygen through different tissues, e.g. the abdominal wall of the frog, or the number and distribution of capillaries in the horse gastrocnemius, and Krogh studied the contraction and opening of capillaries in resting and stimulated muscle [132–134].
As the capillaries appear to be hexagonally packed in skeletal muscle (if viewed in 2D light microscopy cross sections), Krogh assumed in  that for each capillary the 2D region of oxygen supply could be described as a circular area, thus resulting in a so-called Krogh cylinder when viewed in 3D. This averaged area was calculated by counting the number of capillaries within the total muscle area using light microscopy on histological sections. For this geometric limitation imposed by the Krogh model, i.e. a circular region with capillary radius r and cylinder radius R and an oxygen partial pressure P(x) at a point x and under the assumption of a steady state (∂P/∂t = 0), where diffusion in longitudinal direction (orthogonal to the capillary cross-section) is neglected (∂2P/∂z2 = 0) and consumption is constant (MP(P) = p), the diffusion equation (A.1) in appendix A can be solved analytically. The relation between oxygen pressure difference in the capillary (P0) and at a certain point with distance x in the circular region P(x) is then described  as with the diffusion rate D. By setting x = R, the maximum pressure difference necessary to supply the whole muscle is calculated.
Krogh's model was a first useful approximation to the regions of oxygen supply of the capillaries, leading to investigations of the influence of relevant anatomical factors for tissue oxygenation, such as tortuosity and capillary distribution [11–14]. However, studies have found that Krogh's cylinder model is insufficient to describe oxygen supply accurately as it results in non-physiological gaps or overlaps and large regions of unperfused tissue when capillary rarefaction is present . Nevertheless, Krogh's model is still used as a simple approach to model oxygen supply to tissue, for instance recently by Jung et al. . Modification of the no-flux boundary condition at the outer border of the Krogh cylinder into an infinite-domain boundary condition (permitting flux with zero net exchange) leads to a decrease in the axial tissue PO2 gradient . Results of such altered model however remain to be validated by experimental data. Thus, Krogh's model can be used for a first assessment of tissue oxygenation, but more complex models need to be considered, as presented in the following.
4.2. Image-based modelling
The advantage of image-based modelling—in contrast with modelling using simplified and/or artificial geometries—is that the model input represents the actual microvasculature in a realistic fashion. For microvascular networks that have a complex 3D structure, it is difficult to get valid predictions from simple artificial networks. This is even more pronounced in pathological cases such as tumour networks.
To achieve a more physiological description of capillary oxygen supply in comparison to the Krogh cylinder, the concept of ‘domains of influence’ has been introduced. Here, the region of interest is tessellated, i.e. by using Voronoi polygons, which account for heterogeneity of the capillary spacing [137,138]. A Voronoi polygon Vi is defined as the set of points in an area/volume that are closest to a capillary with centre at point xi and farther from any other capillary in the area/volume. Figure 6 shows an example of the Voronoi tessellation (black lines) of a mouse soleus muscle. For tissue tessellated by such supply domains Vi with areas Ai, Hoofd  has generalized the Krogh cylinder solution for the oxygen diffusion to determine tissue oxygen partial pressure at location r to with ri the location of the ith capillary and rci the non-dimensionalization distance, e.g. the capillary radius . Ф(r) is a so-called background function, which depends on the geometry of the oxygen demand domain or domains of influence, such as the Voronoi polygons, and satisfies Green's theorem. Intuitively, the description of the capillary oxygenation regions by Voronoi polygons seems straightforward; under the assumption of an isotropic/homogeneous muscle (diffusion and consumption equal everywhere) and equal and constant PO2 in all capillaries. At all times the oxygen would diffuse across the muscle at the same speed everywhere, thus supplying the regions closest to each capillary and no other, and thus describing Voronoi polygons. Al-Shammari et al.  showed that the description of supply domains using Voronoi polygons is indeed sufficient in most cases. An important exception are pathological cases, where the number of capillaries is greatly reduced (capillary rarefaction), which are not packed quasi-hexagonally within the muscle. By solving the diffusion equation for the 2D steady-state case with constant consumption, Al-Shammari et al.  proved that the areas enclosed by the Voronoi polygons matched the so-called trapping regions well. The trapping regions are defined as those regions where there is no oxygen partial pressure difference across the region boundaries, i.e. the regions supplied with oxygen by a certain capillary, delimiting the streamlines initiating from each capillary. Figure 6 shows the comparison between streamlines (red) and Voronoi polygons (black). The streamlines were generated from the numerical solution to the image-based model (in 2D). The centres of the Voronoi polygons were defined as the geometric centres of the capillaries. In the perfect case, where the Voronoi polygons match the trapping regions (Al-Shammari's hypothesis), no streamlines would cross over the lines of the Voronoi polygons. Al-Shammari's hypothesis was further supported by examining the difference between Voronoi and trapping regions when taking fibre type, size, distribution and varying oxygen demand into account . It was found that the difference in Voronoi and trapping regions was small in most pathological cases, thus making Voronoi tessellations the best and most easily computable representation of capillary oxygen supply regions. However, an investigation and verification of the usefulness of Voronoi tessellations to describe oxygenation volumes in 3D problems has not yet been conducted and remains to be shown.
Most image-based models of tissue oxygen perfusion are based on 2D image data, and namely on histological images [65,139,142]. This is mainly due to the lack of relevant 3D datasets and/or efforts to reduce computing time. 2D models generally come along with the assumption that oxygen diffusion in the longitudinal direction (z-direction) of the muscle can be neglected, to simplify the diffusion process into a 2D problem [139,140,143]. This is a valid assumption for straight vessels along the muscle's main axis in z-direction with parallel cross sections in the xy-plane of the tissue, as the drop of capillary oxygen partial pressure is very low compared with the gradient perpendicular to the vessel. A 3D description can then be obtained by stacking the 2D solutions for all values of z . However, if tortuosity of the vessel is high, i.e. the vessel is sometimes perpendicular to the muscle's main axis, diffusion in the z-direction takes place and can no longer be neglected. An example for this situation was analysed by Penta & Ambrosi . They created a 3D multi-scale model, i.e. capillary intersection (40 µm) versus whole network (1 cm), for drug delivery in tortuous tumour networks, to investigate the influence of capillary tortuosity for drug delivery. Capillary tortuosity was introduced by a sinusoid wave making up the centre line of the capillary, while frequency and amplitude of the wave were varied to increase tortuosity. Hydraulic conductivity and surface area were determined for a microscale cell and then fed into the macroscale problem, which has been solved assuming local periodicity of the microscale cell. In reality, this assumption may not apply, especially for irregular networks like tumour networks. However, the research performed by Penta and Ambrosi shows that tortuosity has a significant impact on drug concentration within the tissue. They showed that the drug concentration in the tortuous network was decreased by more than 50% in comparison to the straight network and that it was below 10% of the injected drug concentration, thus leaving only a small amount to reach the tumour. While only an example case was presented in , this research underlines the importance of treating a vascular network as a 3D structure, and the necessity to take into account the resulting complexity of molecular convection and diffusion. In  the studied network was artificially created and in order to obtain realistic predictions of concentration distributions (here of the drug; in our case oxygen) the computation would ideally be performed on experimentally derived images of tumour vasculature, and the results would need to be validated against experimental data. In a novel approach, Grimes et al.  investigated tumour oxygenation by deriving an oxygenation kernel model, allowing straight vessel segments to diffuse oxygen in a spherical fashion as an array of discrete point sources. Two spherical, scalable kernels were computed; one for a well-perfused tumour network and one for a poorly perfused one, i.e. best and worst case scenarios. The model is inspired by other similar problems, e.g. traffic flow, and is applied to 3D images of tumour vasculature obtained using two-photon microscopy. Validation using fluorescence microscopy was attempted, staining for hypoxia in tumours of the same mouse model. The relative areas of hypoxia resulting from the two kernel models covered a broad range: 37.2 ± 9.8% for high vessel perfusion versus 53.5 ± 13.2% for low perfusion. While the experimentally determined area was in that range (39.6 ± 9.1%) it remains inconclusive how well the kernel model approximates tissue oxygenation as the vessel perfusion level is unknown. The authors argue that validation using 2D methods may be insufficient, but it is also the case that the z-resolution of the applied 3D imaging method was too low to capture all blood vessels and therefore tissue oxygenation may be underestimated. It will be interesting to test the presented approach on a more regular vessel network (e.g. in healthy muscle) captured at higher spatial resolution and using a 3D imaging method, such as light sheet fluorescence microscopy, for validation. Furthermore, while some initial analysis has been provided by Grimes et al., additional analysis needs to be performed to examine the difference between the kernel method and a full computation of tissue oxygenation using a finite-element model.
The works of Goldman and colleagues [12,13,146–148] are based on histological images of the capillaries in hamster cheek pouch retractor muscle, which was subsequently adapted to rat EDL muscle. To obtain a 3D network the 2D information was extrapolated into 3D. The effect of anastomoses and tortuosity was investigated by overlaying straight tissue segments with a sinusoidal wave and adding a physiological number of anastomoses in a random fashion . The authors argue that tortuosity with the addition of anastomoses significantly increased tissue oxygen heterogeneity; however, no tests of statistical significance were provided, thus leaving the results inconclusive. Goldman et al. investigated tissue oxygenation during sepsis [146,147] and it was predicted that pathological conditions were influenced by heterogeneous shunting of blood flow between capillaries. More recently, the same group also compared the oxygen delivery of a 3D reconstructed microvascular volume from intra-vital video microscopy against that of parallel capillary arrays and found significant differences in resulting tissue PO2 . In the healthy and resting case the parallel network has overestimated the PO2 by around 7%, yet, this starkly increased to 18% in hypoxia and 37% during exercise. This demonstrates the importance of realistic 3D microvascular networks over artificially created ones.
Overall, while the mathematics of oxygen diffusion in skeletal muscle tissue are well established and frequently used as in the models presented above, to the best of our knowledge, no image-based studies exist where conclusions have been drawn about links between modelling results and structural parameters of the imaged vasculature. Furthermore, most models include two major limitations: neglecting the influence of the haematocrit and lacking validation against experimental data.
4.3. Limitation 1: role of the haematocrit
Many 2D image-based models of oxygen delivery in skeletal models do not take into account the haematocrit distribution, since they are lacking the structural information of the network of the blood vessels present in the experimentally derived image datasets. Instead, a fixed tube haematocrit is assumed and the discharge haematocrit computed using an empirically derived formula . In addition, the blood oxygen level is kept constant , which is another constraint. Furthermore, it is common to neglect the fact that not all capillaries contain RBCs and thus do not deliver oxygen to the tissue. In other words, the presence of a capillary is equated to the presence of an oxygen-supplying RBC , although the presence of a capillary is only a necessary but not a sufficient condition for oxygen supply. In theoretical mathematical modelling the important role of the RBCs has long been taken into account; some models consider RBCs as point-like oxygen sources , while others assume separate homogeneous flowing regions of plasma and RBCs [149,150]. Gould & Linninger  reviewed such existing models against morphological data of vascular trees and as a result, presented a working model of haematocrit distribution. For image-based modelling of oxygen supply in tissue, this leaves two alternatives: firstly, to image the blood vessel structure and to employ a haematocrit distribution model to correctly predict oxygen delivery; or secondly, in addition to imaging the vascular structure, using an imaging technique capable of identifying RBCs within the vasculature. Tateishi et al.  have imaged the RBCs in an isolated mesenteric blood vessel using an inverted light microscope. They were able to simultaneously measure oxygen saturation and distribution of RBCs within the blood vessel. However, this approach was limited to 2D and to an isolated environment. For 3D imaging in vivo of RBCs, Kamoun et al.  proposed the use of CLSM or multi-photon laser scanning microscopy (a method similar to CLSM, however not using a pinhole aperture but the laser beam itself for optical sectioning of the sample) to capture fluorescently labelled RBCs. Phase-contrast SR CT can also be extended to image the RBC distribution within a 3D tissue sample; however, the information will then still be intrinsically ex vivo, but closer to the in vivo case than in any existing image-based tissue oxygenation model. We have recently trialled imaging of RBCs in mouse soleus muscle using SR phase-contrast CT. Slices of the reconstructed volume are shown in figure 7. For a small cubic region of interest the tissue oxygenation in resting mouse muscle, based on the present RBCs, was computed, see figure 8. Figure 9 shows an excerpt of the 3D visualization of the segmented RBCs and bigger blood vessels.
4.4. Limitation 2: validation of results
The validation of computational modelling results is often not properly addressed. The only mathematical model presented in this section, validated against experimental data, is the tumour oxygenation model presented by Grimes et al. . Complexity of experiments and lack of available resources are most likely the reasons for this. The method employed by Grimes et al.  stains for hypoxic regions in tumour tissue; however, for healthy soft tissue this can pose a problem as no hypoxic regions should in theory be present. Instead, a tissue stain to mark oxygenation could be used, allowing validation of the results from computational modelling by use of gold standard light microscopy. For the validation of a blood haemodynamics model Liu et al.  perfused the muscle of an anesthetized rat with Tyrode's solution. Similarly, in an effort to validate computational modelling results for tissue oxygen supply, the vasculature could be perfused with fluorescent nanoparticles and their diffusion into the tissue could be recorded visually, similar to the work of Kamoun et al. . Another, likely more robust approach is to make use of in vivo imaging techniques that are able to capture oxygenation and blood vessel morphology. There are a number of in vivo imaging methods that are already widely used for blood flow and oxygenation imaging, such as laser Doppler flowmetry  or laser speckle contrast imaging . These imaging techniques generally suffer from their usage of arbitrary units, thus allowing comparison between measurements only, but not quantification of blood flow per se. Japee et al.  introduced a video imaging system for capturing RBC flow and oxygenation simultaneously for superficial imaging of excised muscles. It was applied in Ellis et al.  and Fraser et al.  for comparison with computational modelling results. Other researchers combined different imaging systems, such as optical coherence tomography (OCT) and fluorescence and hyperspectral imaging  and two-photon microscopy combined with a PO2 nanoprobe . Both OCT and two-photon (or multi-photon) laser scanning microscopy have been used extensively by the group of Jain for the characterization of tumour vasculature, but are both limited to imaging only a superficial fraction of tissue within a small field of view [3,158,159]. These methods are however well established and could—with addition of an oxygenation-sensing technique such as a PO2 probe—be used to validate results of image-based modelling of muscle oxygenation as images of the vasculature could be cross-correlated. Thus, the combination of existing imaging and oxygen measurement methods can lead to further exciting results to understand the supply of oxygen to tissue in health and disease.
5. Summary and conclusion
Image-based modelling of skeletal muscle oxygenation is a challenging process. At present, only a few models are based on 3D physiological data and proper validation of computational modelling results is largely missing. As the microvasculature is a complex 3D network its structure needs to be assessed in appropriate detail to obtain information about its capability to deliver oxygen. A wide range of imaging techniques for soft tissue imaging exists, but only few provide the capability to capture small capillaries at whole-organ scale. Structural parameters that are currently in use to investigate this capability in health and disease need to be assessed by mathematical modelling of tissue oxygenation to evaluate their relevance.
Following a short introduction to the biology of oxygen delivery in skeletal muscle, we have reviewed the most common of these structural measures used to assess muscle oxygenation. We have highlighted the suitability of synchrotron radiation-based computed tomography and laboratory-based micro-computed tomography to do this and reviewed necessary methods to achieve image contrast in soft tissue. We further discussed mathematical models of muscle oxygenation with particular emphasis on image-based models. The major limitations of such models, in conjunction with the ex vivo characteristic of the underlying imaging method, have been identified. It becomes clear that further work needs to be done in order to better incorporate the effect of haematocrit and varying capillary PO2 in mathematical models, as well as in finding methods to validate theoretical findings. We have therefore highlighted the synergies that can be gained when combining in vivo imaging methods with oxygen-sensing methods to validate mathematical models and to further progress the assessment of microvascular oxygen exchange capability. Evidently, the emergence of imaging systems at higher spatial resolutions is a great advantage in understanding the link between microvascular structure and its oxygenation potential in health and disease.
B.Z.P. drafted the manuscript and created all original images. T.R., G.F.C. and P.S. helped to draft, and revised, the manuscript. All authors gave approval for final publication.
We declare we have no competing interests.
Funding was provided by the Engineering and Physical Sciences Research Council (EPSRC) for the EPSRC doctoral training grant of author B.Z.-P. and the British Heart Foundation for grant no. PG/12/18/29453.
The preliminary data of red blood cell imaging using phase-contrast synchrotron-based computed tomography were obtained at the TOMCAT beamline of the Swiss Light Source (SLS). We acknowledge the Paul Scherrer Institut, Villigen, Switzerland for provision of synchrotron radiation beamtime at the TOMCAT beamline of the SLS and would like to thank Pablo Villanueva and Alessandra Patera for their assistance. Finally, we would like to thank Prof. Ian Sinclair for valuable discussions that enabled our work and the µ-VIS X-ray Imaging Centre at the University of Southampton for providing the tools for image processing and meshing.
Appendix A. Modelling oxygen diffusion
Oxygen delivery to tissue is driven mainly by diffusion, which is the process of molecules filling a space through random movement. In the steady state this process can be described relating flux to concentration. The relation of flux and concentration is described by Fick's first law with C the oxygen concentration, D its diffusion coefficient (in muscle tissue) and J the oxygen flux . The concentration C(x, t) is dependent on space and time t. The diffusion coefficient D(T) is in fact a matrix, which may be dependent on the direction and the temperature T of the muscle tissue. For simplicity, it is usually assumed that the diffusion process is isotropic and constant and thus the matrix D is a scalar. The diffusion coefficient is normally given at a normal body temperature of 37°C and may vary between different intramuscular soft tissues, such as muscle fibre and interstitial tissue. It is commonly assumed that this variation is negligible and a constant value can be used instead.
The relation of concentration C and partial pressure P can then be described using Henry's law with α being the coefficient of water oxygen solubility . It is assumed in this case that both blood plasma and muscle tissue consist mainly of water and therefore the approximation of α as the water oxygen solubility coefficient is sufficiently accurate.
The change in oxygen concentration at any point in the muscle tissue over time is given by the spatial change in oxygen flux and the oxygen consumption at this point in space. This is the process described in the schematic in figure 3.
The conservation of oxygen is then given by the diffusion-reaction equation or in terms of the oxygen's partial pressure  A 1where M and MP are the respective oxygen consumption functions of the skeletal muscle tissue. They may depend on the oxygen concentration/partial pressure.
In order to solve the above equations it is necessary to employ initial and boundary conditions. As the time derivative ∂P/∂t is first order, only one initial condition is required. On the other hand, is of second order and requires boundary conditions at all muscle tissue boundaries, i.e. here at the vessel–tissue boundary and at the outer tissue boundary.
A.1. Boundary and initial conditions
The flux of oxygen from the capillaries into the tissue is described by a flow along a concentration or partial pressure gradient [12,13,139], depending on the vessel permeability (process (2) in figure 3). The permeability of the wall to oxygen molecules can be described by a mass transfer coefficient k, which depends on the blood haematocrit, the blood flow rate and the type of vessel, i.e. capillary, arteriole or venule, and their radius . This boundary condition is given by where nv is the unit normal vector to the vessel boundary, P and P0 are the oxygen partial pressures at time t and at the onset (t = 0), respectively, and D is the diffusion coefficient. In most cases, k is assumed constant and given by an experimentally determined value (from the literature). This approximation will be valid in most cases; however, in pathological cases where endothelial function is restricted k may play a major role for oxygen supply and must thus be determined. A number of mathematical models have been employed to investigate the influence of different factors on the permeability coefficient, but no uniform formula has been proposed to date. For more information see the review of Goldman .
The boundary conditions are usually imposed at the outer tissue boundary as a no-flux boundary condition: where n0 is the unit normal vector to the outer tissue surface. This means that all oxygen supplied by the blood vessels in one muscle stays within this muscle. This will usually be the case, as diffusion distances for oxygen are small.
The initial conditions for the model can then be set to with Prest being the muscle partial oxygen pressure at rest.
A.2. Oxygen tissue metabolism
The consumption of oxygen within tissue ((3) in figure 3) is usually described using Michaelis–Menten kinetics: where M0 is the maximum oxygen consumption rate within the tissue and P50 the oxygen pressure at half-maximal consumption. P50 is usually found to be 0.5–1.0 mmHg  and can be determined experimentally by methods described in . For skeletal muscle, Al-Shammari and co-workers  took into account the fact that M0 can differ between different muscle fibre types (I, IIa and IIb). The use of Michaelis–Menten kinetics is only an approximation of the occurring process, describing it more accurately than any other model .
A.3. Myoglobin-facilitated oxygen transport
The myoglobin (Mb) oxygen saturation is determined using the Hill equation : with PMb,50 being the half-saturation value of myoglobin. The shape of the myoglobin oxygen saturation curve and of the corresponding haemoglobin (Hb) saturation curve is shown in figure 10. For the same partial oxygen pressure the haemoglobin O2 saturation is lower than that of myoglobin, which indicates that haemoglobin is more likely to release oxygen. Myoglobin-facilitated transport reaches a plateau for a tissue oxygen partial pressure higher than PO2 = 2.5 mmHg and can thus be neglected .
When considering myoglobin-facilitated transport of oxygen ((4) in figure 3), the right-hand side of the general diffusion equation is extended by the term , with DMb, CMb and SMb being the myoglobin diffusivity, oxygen binding capacity and oxygen saturation, respectively. The diffusion equation is thus transformed into
A.4. Haemoglobin oxygen saturation
The oxygen saturation of haemoglobin (SO2) is computed using the Hill equation [130,131]: with the PO2 at 50% haemoglobin saturation, which is determined experimentally, e.g. by using PO2 electrodes on blood samples  or by microspectrophotometry . Together with the haematocrit distribution, this needs to be taken into account for modelling tissue oxygenation, because not all RBCs may be fully saturated with oxygen. n is obtained experimentally and has been reported to be 2.2 for hamsters and 2.7 for rats . Using the Hill equation for modelling haemoglobin oxygen saturation is an approximation of the inverted Adair equation, which is realistically describing the oxygen–haemoglobin dissociation curve .
A.5. Model parameters
The parameters necessary for modelling soft tissue oxygenation usually used for skeletal muscle are shown in table 3. The presented parameters were taken from the sources indicated and were determined experimentally (table 3).
- Received December 8, 2016.
- Accepted January 25, 2017.
- © 2017 The Author(s)
Published by the Royal Society. All rights reserved.