Support for the model’s physiological assumptions

VEGF is a factor with a long history. First identified as a vascular permeability factor [30, 31], it has become clear over time that it has other important roles as a factor in angiogenesis [32], endothelial cell proliferation [33], and also as a neuroretinal protective factor [34, 35]. While we do model retinal edema [36] the important property of VEGF for this model is its role a mediator of elevated ICAM-1 on retinal endothelial cells. In our model of progressive capillary occlusion, VEGF is the locally secreted molecule which diffuses and increases the likelihood of nearby capillary occlusion. As this pro-occlusive property is not as well-known and even possibly denied by some, we provide below the background supporting our choice of VEGF as a substance responsible for the adverse cycle of capillary occlusion. This does not mean VEGF is by any means the only cytokine involved or that there is not an intervening cascade of events that generate the chronic inflammatory state that is diabetic retinopathy. Leukostasis is mediated by the diabetic activation of circulating leukocytes co-existing with marked upregulation of adhesion molecules such as ICAM-1 on the retinal vascular endothelium [37]. These changes increase the likelihood of leukocyte adhesion to the retinal capillary endothelium and therefore the probability of capillary closure. The model does not detail the local, undoubtedly complex, phenomena such as the cumulative leukocyte mediated endothelial capillary damage resulting in endothelial replicative exhaustion culminating in capillary occlusion; they are treated as black boxes at this time.

The model treats elevation of one substance for simplicity though likely the relative amounts of two substances, such as the balance of VEGF and pigment epithelium derived factor (PEDF), is what is often physiologically important [38]. PEDF itself has complex neurotrophic, neuroprotective, and anti-angiogenic, anti-exudative and anti-inflammatory properties [39, 40]. High glucose decreases expression of PEDF in retinal Mueller cells as it simultaneously elevates VEGF expression. Additionally vitreous levels of PEDF are significantly lower in patients with diabetic macular edema or proliferative diabetic retinopathy than in non-diabetic patients or diabetic patients without retinopathy whereas in each situation VEGF is elevated. The model simplifies this duality by treating the physiological import of an imbalance as simply the concentration of VEGF. Certain steps in the progressive ischemic process must be met by a diffusible substance and ideally we need to have anatomical and physiological support for each of the steps in the model if we are to model the substance as VEGF. The local adverse positive feedback model could stand on its own dealing only with the problem of diabetic ischemia as a geographic phenomenon but it is more constructive to have model elements that correspond as closely as possibly to biological elements. We now construct a simple biological model based on the observation that this substance is required to have a number of physiological properties which first create the permissive diabetic state of recurrent capillary endothelial cell loss by activated leukocytes, leads to permanent capillary occlusions from local endothelial depletion, and which in turn produces geographic propagation of capillary occlusions. This progression requires a substance, modelled as VEGF, and the retinal tissue to have the following characteristics:

1. Step 1: With the onset of hyperglycemia or the continued presence of hyperglycemia this substance must change its production in some way from that of the non-diabetic state.
2. Step 2: The change in concentration of this substance must cause some change in the retinal capillaries resulting in an increased probability of permanent capillary occlusion above essentially zero in the non-diabetic state.
3. Step 3: Capillary occlusion must in turn create local retinal ischemia and by some mechanism further elevate the level of the substance, meaning that ischemia cannot act to kill the cells producing the substance.
4. Step 4: Higher concentrations of this substance in a local retinal area must cause a higher likelihood of occlusion in nearby capillaries with the probability of occlusion functionally related to the concentration of the substance.

Human subjects

The protocol for AOSLO imaging used to provide retinal vascular images used in this study was approved by the Indiana University Institutional Review Board and adheres to the tenets set forth in the Declaration of Helsinki and the Health Insurance Portability and Accountability Act regulations. Written informed consent was obtained from all subjects.

Model considerations

To capture the events of progressive capillary occlusion we implemented a quantitative model of the anatomical features mentioned above in Compucell3D [98]. As shown in Fig 2, beginning with an AOSLO scan of the perivascular fovea, a full vascular model from an arteriole to a venule with the linking capillaries is included with oxygen advection, oxygen diffusion, and oxygen consumption. This paper develops a conceptually simple model of the diabetic retina treating Mueller cells as the sole retinal source of VEGF and assumes a slight elevation of VEGF production by Mueller cells in a diabetic retina higher than that in the normal retina [99]. Physiologically VEGF is a necessary neurotrophic factor in the retina and is normally present at low levels [100–102]. In the model VEGF is produced by Mueller cells locally in variable amounts based on oxygen saturation. VEGF diffuses from the Mueller cells and is consumed by cells including endothelial cells but is not transported away by advection. The model vessels have endothelial cells which respond to local VEGF levels by an increased probability of occlusion with elevation of local VEGF and also by leaking if local VEGF exceeds a threshold level. The model is cycled many times and if a capillary occlusion occurs, all flow rates, steady state oxygen tension and VEGF levels are recalculated. The model’s treatment of occlusion is an irreversible decrease of capillary diameter to zero. An important assumption is that the vascular supply to each area of retina is critical in that occlusion of a capillary will result in ischemia of an area of physiologically dependent retina with a resultant elevation of VEGF synthesis by the locally ischemic Mueller cells. We do not know that this has been proven but the constraints imposed by evolution on the visual apparatus make this assumption reasonable. To quote Chan et al. 2012 [102] “It is likely that retinal capillary networks are morphometrically adapted in order that the balance between cellular nutrition and optical clarity can be achieved.” Note, however, that with the variations of capillary spacing seen anatomically, all areas of retina would not have equal dependence on a single supplying capillary. In a network based on actual capillary anatomy different areas of retina could be more or less critical as a result of variation in local capillary density. There would thus be greater or lesser propagation of capillary closure by the adverse feedback mechanism. Small capillary diameter adjustments can also occur, e.g. slightly increasing diameter with increased flow after each capillary occlusion. Maps at various model times are made of capillary network structure, flow, oxygen tensions, VEGF, and retinal edema. These are the output measures as well as summary graphs of the system such as total flow and average distance from an intact capillary.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Model abstraction and construction of human perifoveal capillary network from an Adaptive Optics Scanning Laser Ophthalmology (AOSLO) image.

On the right, AOSLO image shows juxtafoveal capillaries adjacent to the foveal avascular zone (FAZ). This is a normal capillary map in a patient without diabetes. On the left, a model schematic shows the reconstruction of the capillary network framed in the AOSLO image with cells filling in empty space between vessels uniformly. Capillary network has arteriole (A) and venule (V) termini marked. Boundary blood pressures are assigned for A and V termini. Boundary oxygen tensions are assigned for A terminus and for FAZ whereas venous oxygen tension is model dependent. Model objects in red, green, brown are capillary blocks (CAP), Mueller cells (MC) and other retinal cells (OT) respectively. Yellow pixels surrounding objects are object borders, which are muted in other figures of the manuscript. Note that this example shows a macular capillary network from a different subject than that shown in CASE 1. Diameter of the FAZ in this ASOLO image is approximately 500 microns. The scale bar is 100 microns.

https://doi.org/10.1371/journal.pcbi.1004932.g002

Model capillary occlusions are always probabilistic based on local VEGF levels and the calculated flows of the capillary segments. Capillary networks of several types were utilized including physiologically unlikely hexagonal capillary network with introduced deletions, physiologically realistic peripheral retinal ‘ladder’ capillaries [103] and an actual perifoveal arteriovenous sector capillary map obtained from adaptive optics scanning laser ophthalmoscopy (AOSLO) imaging of a subject. The hexagonal map was used to explore the dependence of capillary occlusion progression on amount of tissue dependent on capillaries by varying the scale of the hexagons. Both the macula and the peripheral retina are clinically important, with the macular area being the location of ischemia as well as macular edema affecting visual acuity, and the periphery being the major source of the ischemia and resultant VEGF production which results in retinal neovascularization. Though it is not clear if it is physiologically appropriate to utilize the same fundamental model parameters for the macula and periphery given the distinct peripheral retinal and vascular architecture [103, 104] we altered only the capillary diameters and network structures. In the main text we address primarily a sector of the perifoveal capillary network from AOSLO imaging (CASE 1) and filled the open space between vessel segments with cells of anatomically reasonable sizes (S4 Table). Capillary diameters were estimated based on the AOSLO image. Model inputs such as terminal hydrostatic pressure and arteriolar blood oxygen tension were estimated from published results (S4 Table, S15 Fig). Vascular flows, oxygen and VEGF fluxes were calculated and resulting tissue oxygen tension and VEGF levels were determined in the model. A large number (362) of replicate runs were made of the subject’s capillary network in order to assess the vulnerability of distinct capillaries given the probabilistic nature of the model of individual capillary occlusion. S1 Text details mathematical descriptions, parameter selection and influence and boundary and initial conditions of the model. S2 Text briefly treats a run of the model in the macular area with a different initial capillary occlusion (CASE 2). S3 Text treats a hexagonal network with some deleted edges, and S4 Text treats a peripheral schematic ladder capillary network modelled on that of the peripheral human retina. The peripheral ladder capillary network seen in the peripheral retina has significant implications about the development of the peripheral ischemic wedge-like sectors seen clinically between arterioles and venules whereas the hexagonal model was informative about the dependence of capillary occlusion propagation on distance between capillaries.

Many models have been constructed to study problems at the interface of vasculature in various tissues: skeletal muscle [105, 106], brain [107, 108], vascular tumor [109, 110] and retina [111]. Shirinifard et al. employed a 3D multi-cell model to successfully recapitulate the three patterns of progression of age-related macular degeneration and suggested that defects in adhesion were the dominant contributor to initiation and development of choroidal neovascularization [111]. Cringle et al. divided retina into multiple layers and used a mathematical model to calculate the oxygen tension in each layer in terms of oxygen consumption rate in that layer and the oxygen level in choroidal capillaries [112, 113]. McDougall et al. studied angiogenesis during normal retinal development using a hybrid discrete-continuum mathematical model and computationally simulated the structure of a retinal vascular plexus that agreed with the whole-mount retinal vasculatures at different stages of development [114]. Our model deals with a different pathophysiological issue: progression of ischemia and edema in diabetic retinopathy based on a local VEGF-dependent mechanism of propagation of capillary occlusions. Unlike this study, Gandica et al. [115] developed a computational model of retinal ischemia studying the effect of critical sizes and densities of localized blockages of retinal vasculature on the emergence of diabetic retinopathy. In their model, various sizes of local blockages of vessels, assumedly caused by destabilizing proteins such as Angiopoietin-2, were randomly distributed in the region of interest and areas of derived hypoxia were examined regarded as an indicator of potential phenotypes of diabetic retinopathy. An important conclusion in their study is that local blockages with smaller size than characteristic irrigation length, if their densities exceed a critical threshold, likely result in large hypoxic areas because of a cooperating effect. A limitation of the model, as the authors also noted, is the simplified consideration of oxygen transport.

To ensure that we respect the specific geometry of the retinal capillary network, we use networks from actual subject imaging or from peripheral retinal capillary networks in the literature and were therefore not created for these modelling purposes. Also the oxygen diffusion coefficient was taken from the literature.

Model description

In this study, we explore the effect of focal capillary occlusion on decreasing local oxygenation of retinal cells, resultant elevation in VEGF, and the consequences in terms of propagation of capillary occlusions and formation of edema using a computational model. The anatomy of this capillary network was determined from a normal patient using AOSLO (Fig 2). This network is an arteriovenous sector with the capillaries connecting the arteriole and venule, the foveal avascular zone on one edge and retinal tissue on the other borders. We computationally reconstruct the network and initialize the simulation with boundary blood pressures and arterial and FAZ oxygen tension (Fig 2). The major entering arterial node and the exiting venous node are assigned with blood pressures and respectively. These pressures would not be expected to depend on the diabetic state and remain fixed. Other boundary nodes are assigned with blood pressures of intermediate values (S15 Fig, also see Boundary conditions and initial state of simulation in S1 Text for details). The FAZ region is treated as an oxygen source, supposedly supplied by choroidal capillaries. The entering arterial node and the FAZ region are assigned with oxygen tension and though the exiting venous oxygen tension is model dependent. Other boundary nodes with incoming blood flow are assigned with oxygen tensions of smaller values (S15 Fig). Hypoxia-induced elevated secretion of VEGF in Mueller cells is known as an occurrence in diabetic retinopathy [115–117]. FAZ region is treated as a sink for VEGF. We modelled the oxygenation and local VEGF levels within this retinal sector following a focal capillary obstruction. Within the sector of capillary network modeled, the arterial side receives oxygen-rich blood, while from the venous terminus carries away blood with lower oxygen tension. Oxygen diffuses into tissue space from capillaries, where it is consumed and metabolized by retinal cells. Hypoxia of local retinal tissue is induced by a local capillary segment occlusion, and production of local VEGF, dependent on the level of hypoxia is upregulated in ischemic Mueller cells. VEGF released by Mueller cells diffuses to other nearby patent capillary segments, which probabilistically derive more occlusions based on an underlying and not visibly modelled upregulation of ICAM and increased probability of leukostasis. The schematic of oxygen and VEGF fluxes is shown in Fig 3.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Schematic of oxygen and VEGF fluxes.

Colored blocks represent model objects: Capillary block (CAP in red), Fluid portion (FP in cyan), Mueller cell (MC in green) and Other retinal cells (OT in brown). Markers in the form arrows, triangles and squares are used to represent modeled fluxes. On the left are the oxygen fluxes including advection, diffusion and metabolism. Oxygen advection is modeled for object CB. Oxygen diffusion is modeled for object pairs among CAP, MC, OT and FP. Oxygen metabolism is modeled only for objects MC and OT. On the right are the VEGF fluxes including synthesis, diffusion and decay. VEGF synthesis is modeled for object MC again as the model’s only VEGF source. VEGF diffusion is modeled for object pairs among CAP, MC, OT and FP. VEGF decay is modeled for objects CAP, MC, OT and FP. The arrangement of model objects do not necessarily reflect detailed configurations constructed in the simulations.

https://doi.org/10.1371/journal.pcbi.1004932.g003

Model objects and processes

The described computational model consists of four generalized model cell types: capillary block (CAP), fluid portion (FP), Mueller cell (MC) and other retinal cell (OT). Along with these four generalized cell types in the model is another model object called the conveyor-belt block, CB, which is an object associated with the capillary block and introduced for modeling of oxygen advection. In addition to five model objects, two chemical fields exist in the model: oxygen and VEGF. Modeled processes include advection of blood carrying oxygen, diffusion and metabolism of oxygen, and synthesis, diffusion and decay of VEGF. As detailed in S1 Text and summarized in S1 Table, model objects have the following properties and are representative of various retinal cells:

1. CAP is the structural element of the capillary segment. CAP functionally represents a capillary segment with endothelium and blood within the capillary lumen. CAP is involved in modules of oxygen diffusion and VEGF diffusion and decay. CAP has three states: “normal”, “leaky” and “occluded”. “Normal” CAP can transition to “leaky” state if local VEGF level is greater than a threshold value. Both “normal” and “leaky” CAPs can transition to an “occluded” state if the occlusion condition, whose probability is raised by an elevated VEGF level and lowered by elevated blood flow velocity, is met (S2 Table). Once one CAP becomes “occluded”, the flow of that capillary segment from node to node is set to zero. Also, “Occluded” is a permanent or irreversible state for the CAP. (More details on capillary occlusion in S1 Text)
2. FP is the structural element of and represents incremental volume of retinal edema, cystic areas of fluid within the retina. FP represents fluid leaked from a CAP while it is in the leaky state. FP is modeled as a “cell”-like object that is well trapped by surrounding objects. FP is involved in modules of oxygen diffusion and VEGF diffusion and decay. (More details on edema formation in S1 Text)
3. MC is involved in modules of oxygen diffusion and metabolism and VEGF synthesis, diffusion and decay. The MC has two states: “normal” and “hypoxic”. The “Normal” MC can reversibly transition to a “hypoxic” state synthesizing VEGF only when cellular oxygen tension is below a threshold value. In the model only the “hypoxic” MC has the capability of synthesizing VEGF though in the diabetic retina a number of other cells have at least some ability to synthesize VEGF.
4. OT broadly includes all other retinal cells other than Mueller cells, i.e. the neural retina including astrocytes and microglia but excluding photoreceptors supplied with oxygen from the choriocapillaris. OT is involved in modules of oxygen diffusion and metabolism and VEGF diffusion and decay.
5. CB is a functional element created in the model to “convey” oxygen along a linear pipe, thus CB or conveyor belt is not an anatomical element corresponding to a retinal structure but is a model device allowing mathematically convenient modelling of oxygen advection. CB is involved in the module of oxygen advection. Each CB is associated with a “host” CAP, which usually contains more than one CB. In contrast to the fixed pre-defined size of CAP on every capillary segment, size of CB is proportional to the blood flow velocity on its host capillary segment.

To model oxygen advection, we discretize each capillary segment into a one-dimensional sequence of equally-sized CBs and simulate oxygen advection using a “conveying” action, which moves a volume of oxygen in a given CB to its next downstream connected CB (Fig 4). The size of a CB is proportional to both flow velocity and the time step for advection. In the case of merging at a junction, total oxygen volumes in last CBs of the predecessor capillary segment are conveyed to the first CB of the successor capillary segment (Fig 4). In the case of branching at a junction, conservation of blood flow is enforced to appropriately distribute the oxygen volume in the last CB of the parent capillary segment to the first CB of the daughter capillary segments (Fig 4). Mathematical descriptions were detailed as equations (9)-(11) in S1 Text. The benefit of the conveyor-belt model of oxygen advection is the flexibility of adjusting the number of CBs in a CAP without interference with oxygen diffusion. When occlusion occurs, flow velocities change on the capillary network and accordingly the size and number of CB on a CAP also change. This alters the discretization for oxygen advection, but the module of oxygen diffusion remains unaffected. A simple rule is followed to associate a CB with CAP: a CB is associated with a CAP if the CB’s center is enclosed by the CAP’s volume. Then, the following sequence of processes is used to link modules of advection and diffusion: (1) oxygen advection in CBs on each capillary segment at time step t₀; (2) conversion of oxygen volume from CB-level to CAP-level just before diffusion at time step t₀; (3) oxygen diffusion involving CAPs and their surrounding objects at time step t₀; (4) conversion of oxygen volume from CAP-level to CB-level just before advection at time step (t₀ + Δt_f) (Fig 5). (More details on the CB model of oxygen advection are in S1 Text).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Schematic of conveyor-belt model of oxygen advection.

The schematic shows a small capillary network with five segments and two junctions. Along the blood flow direction, marked as arrow with empty triangular head, the first junction has merging blood flow and the second branching blood flow. Blocks in red are CAPs, building blocks, or structural elements (visually present in model configuration), of capillary segments. Blocks in yellow are CBs, functional elements for advection. Size of a CAP is fixed throughout a simulation, so each capillary segment has fixed number of CAPs. In contrast, the size of a CB is always proportional to flow velocity on the capillary segment, so flow velocity and accordingly the number of CBs on a patent capillary segment can vary following an occlusion elsewhere as flows in the network are adjusted to the changed network resistance structure. During each time step of advection, a CB passes its oxygen volume to the next CB. Importantly, the size of a CB is equal to flow velocity of host capillary segment multiplied by time step of advection, which means the “conveying” speed of oxygen from one CB to next is exactly equivalent to the flow velocity in that capillary segment. At a merging junction, the upstream capillary segments add the oxygen volumes in their last CBs and pass the total to the first CB of the down stream capillary segment. At a branching junction, the parent capillary segment distributes the oxygen volume in its last CB into the first CBs of the daughter segments according to conservation of blood flow volume. Mathematical descriptions were detailed as equations (9)-(11) in S1 Text. Each CB is associated with, or mapped to, a host CAP. The modules of advection (involving CBs) and diffusion (involving CAPs) are connected with processes that convert oxygen volumes between the host CAP and its associated CBs. (see Fig 5).

https://doi.org/10.1371/journal.pcbi.1004932.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Procedure of sequential simulation of advection and diffusion illustrated with an example.

During each time step, simulation of advection precedes simulation of other processes including diffusion. After advection (involving the CBs delivering the blood containing oxygen) and before diffusion (involving CAPs), a CAP always sums up oxygen volumes in all its associated CBs and update its pre-diffusion oxygen volume. Simulation of diffusion updates the CAP to have its post-diffusion oxygen volume. After diffusion and before advection at the next time step, the associated CBs must have the same relative change of oxygen volumes as the CAP gains or loses during diffusion, namely, percent change from pre-diffusion to post-diffusion values. This mathematically recognizes the conservation of oxygen. Three consecutive time steps are shown in the example. A naïve assumption here for this illustration is that the first CB always receives 1 unit oxygen volume from upstream (not drawn). Another simplification made for this illustration is that only diffusion between the CAP (shown in red) and nearby OT (shown in light brown) are considered. In the actual model OT to OT diffusion is also treated. Note that in a certain step, boldface numbers represent values being changed or updated. During the first time step from t₀ to t₀ + Δt, while the second CAP and its associated CBs still have a zero oxygen volume, the first CAP and associated CBs undergo (1) the process of advection that passes 1 unit oxygen volume to first CB, while CAPs are not involved; (2) an intermediate step that updates CAP’s pre-diffusion oxygen volume by adding 1 (its first associated CB) and 0 (its second associated CB); (3) the process of diffusion delivers 0.2 to OT in contact (amount assumed for convenience in this example, and again OT-OT diffusion is ignored in this example) and CAP’s post-diffusion oxygen volume becomes 0.8, while CBs are not involved; (4) a last step in the time period that updates CAP’s associated CBs’ oxygen volumes by subtracting 0.2/1 = 20% (diffused/pre-diffusion), the first CB thus having 0.8 oxygen volume now. During the second time step from t₀ + Δt to t₀ + 2Δt, similar verbal “simulation” goes. (1) process of advection goes as another 1 oxygen volume is passed to first CB and 0.8, previously held by the first CB, is passed to the second CB; (2) an intermediate step adds 1 and 0.8 to first CAP, still none added to the second CAP; (3) the process of diffusion updates first CAP’s oxygen volume to 1.44, with 0.36 diffused out; (4) last step updates oxygen volumes in both of the host CAP’s two CBs, again by subtracting diffused fraction 0.36/1.8 = 20%. During the third time step from t₀ + 2Δt to t₀ + 3Δt, advection now passes oxygen volume 0.64, previously held by the second CB, into the third CB, which is associated with the second CAP. An intermediate step updates both CAPs by summing up oxygen volumes in their associated CBs. Process of diffusion now changes the oxygen volumes of both CAPs, with the first and second diffusing out 0.36/1.8 = 20% and 0.1/0.64 = 15.625% respectively. The last step subtracts the oxygen volumes of their asscociated CBs’ with the percent change.

https://doi.org/10.1371/journal.pcbi.1004932.g005

To model diffusion of a chemical field, we assume that each generalized cell has a uniform intra-cellular chemical concentration and that diffusion occurs at the interface between neighboring cell pairs. Metabolism of oxygen and synthesis and decay of VEGF are modeled as an intra-cellular process (Fig 3). As for diffusion, the rate of exchange of oxygen or VEGF between generalized cells is proportional to inter-cellular gradient of concentration and inter-cellular contact surface area. Metabolism of oxygen obeys Michaelis-Menton kinetics, which assumes variability of oxygen consumption given different intracellular oxygen tension. Synthesis of VEGF in hypoxic MC is dependent on intra-cellular oxygen tension and VEGF level. Decay of VEGF is modeled using first-order kinetics. FAZ region is treated as both an oxygen source and a VEGF sink. (Details on oxygen and VEGF fluxes in S1 Text)

Model workflow

Our simulation involves three distinct intervals of time: the time step of integration for oxygen and VEGF flux, Δt_f; the time interval to check edema formation, Δt_e; and the time interval to check for occlusion, Δt_o. The Δt_f is chosen so that the differential equations of fluxes are properly integrated. In contrast, Δt_e and Δt_o are selected so that possible edema formation and capillary occlusion take place at a significantly slower pace (months to years), as compared with fast establishment (seconds) of the oxygen and VEGF steady state following a newly derived occlusion. With Δt_f = 0.002s, Δt_e = 7 days, Δt_o = 28 days, a model workflow is shown in Fig 6. The time for oxygen and VEGF steady state is basically mandated by physics whereas the other longer times for edema and capillary occlusion are chosen to allow the model to progress in rough accordance to the progression of clinical disease.

1. Input of network anatomy. The simulation starts with input of a capillary network and cells. Both structural and topological information of the capillary network is required. In our model, initial structural and topological information is determined by digitization of an experimental image obtained by AOSLO imaging. Diameters of capillaries (d^cap were estimated to follow a Gaussian distribution with mean value of 5 μm and standard deviation of 0.5 μm based on measurement of the AOSLO image, and the terminal venule and arteriole were measured to have larger diameters up to 10 μm. This structural information, with topological information of nodes and edges, is read into Compucell3D to reconstruct the capillary network. Next, the Mueller cells, MCs, were organized uniformly into the open space surrounding the capillary network. These are of roughly anatomical size. The space not occupied by MCs becomes the other tissue (OT). (More details in boundary conditions and initial state of simulation of S1 Text)
2. Computation of flow velocities. Blood flow directions and rates were determined using the Poiseuille equation, with hydrostatic pressures at boundary nodes given and fixed through the simulation. Boundary nodes refer to all those extending to the capillary network outside the region of interest, including the arterial inlet and venous outlet. (More details in Network flow of S1 Text).
3. Simulation of oxygen fluxes. Simulation of oxygen advection, diffusion and consumption are divided into three separate sub-modules sharing the same time step of integration Δt_f. An important feature of the model is that advection is simulated at CB level, while diffusion and consumption are simulated at inter-cellular and cellular levels respectively. Therefore, before actual simulations of advection, a process is needed which reads current CAP-level oxygen volume to update CB-level oxygen volume. After simulation of advection by the CBs and before simulation of diffusion and consumption of oxygen, a process is needed to convert CB-level oxygen volume to CAP-level oxygen volume. The simulation of oxygen fluxes is executed until all model objects arrive at a steady state of oxygen tension. (More details in Oxygen flux of S1 Text)
4. Simulations of VEGF fluxes. Simulations of VEGF synthesis, diffusion and decay are divided into separate sub-modules which share the same time of integration as the oxygen flux, Δt_f. This module follows immediately after the simulation of oxygen fluxes. The simulation of VEGF fluxes is executed until all model objects arrive at the steady state of VEGF level. (More details in VEGF flux of S1 Text).
5. Formation of edema. Every period of Δt_o, the condition for edema formation is checked for all CAP objects. Once local VEGF level at a CAP exceeds a pre-defined threshold , edema is formed near the CAP. The edema is modeled as a pseudo cell FP, which is created every Δt_o as long as the requirement of supra-threshold VEGF is met, namely, the CAP remains in a “leaky” state. Fluid portions of the pseudo cell get eliminated at the bottom of the system, which represents the function of the retinal pigment epithelium removing excess fluid (More details in Edema formation of S1 Text).
6. Probabilistic dependence of capillary occlusion. Every period of Δt_e, the probability function for determining capillary occlusion is checked for all CAP objects. If the requirement, which is dependent on the intracellular VEGF level of a CAP and blood flow velocity of the capillary segment, is met, the CAP turns to the “occluded” state, and so do all CAPs that belong to the same capillary segment. Consequently, the capillary segment acquires an infinite flow resistance (in the model, an infinitesimal number is now assigned to the diameter of the segment) and is effectively occluded. (More details in Capillary occlusion of S1 Text) This is the model’s representation of the process of capillary occlusion by a leukocyte. A high flow is assumed to make adhesion of a leukocyte less likely as mechanically the leukocyte is experiencing more pressure pushing it down the capillary. The dependence on level of VEGF subsumes the dependence of local ICAM levels on local VEGF with higher VEGF levels giving higher ICAM levels and increased leukocyte adhesion.
7. On occurrence or absence of a new capillary occlusion. If occlusion occurs, network topology changes and capillary diameters are slightly adapted in response to hemodynamic and metabolic stimuli. Steps (2)-(6) are repeated until the end of the simulation or until all patent capillary flow paths no longer exists. If no occlusion occurs, steps (5)-(6) are repeated until a new occlusion occurs or until the end of the simulation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Model workflow.

Simulation is initialized with input of the original capillary network topology and structure as well as the component cells. Computation of flow velocities and simulation of oxygen and VEGF fluxes occurs next. Each model week the model addresses capillary leakiness and each model month capillary occlusion is evaluated. (1) If no new occlusion occurs, the model simply repeats checking for edema formation and capillary occlusion until occlusion occurs or until a pre-assigned simulated time (3 years at our temporal conversion rate) is arrived at. (2) Once a new occlusion occurs, the network topology is changed and its structure is adapted to eliminate flow in the irreversibly occluded capillary, followed by a new iteration of previous steps. The process stops if either no more effective flow paths exist, invalidating flow velocity calculations, or the pre-assigned simulated time is arrived at.

https://doi.org/10.1371/journal.pcbi.1004932.g006

Model Outputs

(1) Capillary occlusion.

The model starts with the cell-vessel configuration show in Fig 7A. In order to visually track occlusions of the capillary network, blood flow pathways are shown using a color map at each time point (Fig 7B). Additionally the model calculates the average minimal cell-to-vessel distance d^min, which measures the distance from a Mueller cell to its closest patent, or non-occluded, capillary segment as a system output.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. (CASE1) system configurations and field pattern under normal condition.

(A) Shown in the 2D XY cross section involving capillary network, Mueller cells (green) and other retina cells (brown) are uniformly initialized between capillary segments composed of capillary blocks (red). The 3D image shows the simulated section during initialization (in the 3D completely initialized configuration, the capillary network is covered by MC and OT cells and visually inaccessible). The dimension of the simulated retinal section is 510μm × 600μm × 50μm. (B) Flow velocity map includes a primary arteriolar entrance (bottom right in the (B)), a main venular exit (top right), side traffic extending outside region of interest, and interconnected pathways. Capillaries closer to the FAZ (left in the figure) carry blood flow with relatively smaller velocity. The unit for velocity is μm/s. (C) Oxygen tension is highest near capillary segments and the FAZ. The more distant a cell is from irrigating capillaries, the lower is its oxygen level. The unit for oxygen tension is mmHg. (D) VEGF levels are initially assumed to be at a low but nonzero basal concentration across in the whole area. VEGF level has arbitrary unit. “FAZ refers to foveal avascular zone, “A” refers to arteriole, and “V” refers to venule.

https://doi.org/10.1371/journal.pcbi.1004932.g007

Spatial patterns of disease progression

It is important to maintain perspective regarding the model. While we initially present results in terms of an example run of the model, CASE 1, as explained below, the model was run many times to explore the impact of the implemented stochastic events.

The configuration of cells and vessels has been initialized for a sector near the FAZ with dimensions of 510μm × 600μm × 50μm, as determined from a patient ASOLO image. The sector during initialization is viewed in 3D with CAPs and MCs visualized (Fig 7A). Completely initialized configuration shows that MCs and OTs are uniformly patterned between vessels from a 2D view (Fig 7A). Under normal conditions, the blood flow is sourced from the arteriolar entrance, flows through the capillary network and exits via the venule. Because the model concentrates on a small arteriole-venule sector, side streams reflect capillaries connecting neighboring peri-foveal networks (Fig 7B). Capillaries near the FAZ carried a relatively small blood flow, and therefore had a higher probability of occlusion based on the assumed occlusion mechanism. Cells near capillaries had abundant oxygen supply, while more distantly situated cells received less, but still sufficient, oxygen to support normal activities without ischemia (Fig 7C). In the initial state VEGF production was assumed to be in the basal diabetic physiological range when retinal tissue was adequately oxygenated (Fig 7D). This corresponds to a slight elevation of VEGF above the basal non-diabetic state; an elevation sufficient to induce the presence of a low level of ICAMs allowing diabetes induced leukostasis.

Once a first occlusion was initiated, the network topology changed reflecting the loss of that capillary and the steady state of oxygen tension and VEGF level were reestablished accordingly (Figs 8–10). In week 0, a capillary was occluded (Fig 8A) and consequently tissue surrounding it became poorly oxygenated (Fig 9A). Within the model, the date of first capillary occlusion is always week 0. This is not the date of the initiation of diabetes. Certain Mueller cells become ischemic from this capillary occlusion and these hypoxic Mueller cells produce and release VEGF, which then develops local concentration peaks (Fig 10A). Adjacent capillaries respond to the elevated concentration of VEGF, which increased their risk of occlusion via the ICAM mediated leukostasis mechanism not addressed separately in the model. Due to limited diffusion length of VEGF, distant capillaries are insensitive to such localized change. In week 72, secondary occlusions took place and the flow network became increasingly impaired (Fig 8B). Consequentially, more Mueller cells had an insufficient oxygen supply and produced increased levels of VEGF (Figs 9B and 10B). At this point, the capillary damage appeared well bounded and confined in one arteriolar-venular sector. Notably, several cells in a neighboring area became hypoxic and produced elevated VEGF, possibly because their nearest vessel lost an important upstream branch, which carried significant oxygen supply before closure. However, there were still barriers composed of healthy cells and capillaries between the adjacent sectors. There then followed a series of capillary occlusions with minor enhancement in total volume inflow rate, until the damage in week 124 resulted in drop of inflow rate the first time since onset of initial injury (see below).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. (CASE1) flow velocity pattern following initial capillary occlusion.

(A) The flow velocity map captures loss of a flow pathway due to capillary occlusion in week 0. (B) The second capillary spatially close to initial occlusion site became occluded in week 72. (C) Several capillaries near FAZ and venule became occluded in week 124. (D) More than a quarter of the capillary network was obstructed by week 152. Color and the direction of arrows reflect the magnitude and orientation of velocities respectively. The redder the color is the greater the flow velocity. The unit for velocity is μm/s. “FAZ” in the figure refers to foveal avascular zone, “A” in red refers to arteriole, and “V” refers to venule.

https://doi.org/10.1371/journal.pcbi.1004932.g008

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. (CASE1) oxygen tension pattern following initial capillary occlusion.

(A) Oxygen tension map shows a localized hypoxic region near the occluded capillary in week 0. (B) Hypoxic area of cells broadened to enclose a second occlusion site in week 72, but it’s still restricted and confined spatially to the arteriole-venule sector. (C) Local insult by hypoxia propagated to break the terminal capillary near the venule and extended to the neighboring AV sector in week 124. (D) Large area of hypoxia was observed in week 152. Color reflects magnitude of oxygen tension. The redder the color the higher the oxygen tension is. The unit for oxygen tension is mmHg. “FAZ” in the figure refers to foveal avascular zone, “A” in red refers to arteriole, and “V” refers to venule.

https://doi.org/10.1371/journal.pcbi.1004932.g009

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. (CASE1) VEGF level pattern following initial capillary occlusion.

(A) VEGF level map shows localized synthesis of VEGF by Mueller cells in response to hypoxia in week 0. (B)—(D) Increasing number of Mueller cells actively produced VEGF in weeks 72, 124, and 152, where the pattern of regions with high VEGF reproduced that of the area with low oxygen tension. Color reflects magnitude of VEGF level. The redder the color the higher the VEGF level is. VEGF level has arbitrary unit (A.U.). “FAZ” in the figure refers to foveal avascular zone, “A” in red refers to arteriole, and “V” refers to venule.

https://doi.org/10.1371/journal.pcbi.1004932.g010

Concomitant with the drop in total flow, the terminal venule between the two sectors was compromised (Fig 8C top of image) and a larger fraction of tissue was now hypoxic (Fig 9C). In week 152, capillaries near the FAZ and the terminal venule were no longer patent (Fig 8D), the FAZ is considerably enlarged, and roughly one third of the tissue was hypoxic and had an elevated VEGF environment (Figs 9D and 10D). We note that Sakata et al. 2006 [118] using a fluorescein angiography technique to measure perifoveal capillary blood velocity were able to show a significant negative correlation between capillary blood flow velocity and the size of the foveal avascular zone in diabetics without edema. In the model, capillaries bordering the avascular zone show slowed flow consistent with Sakata et al. results [118].

Edema formation

Based on simple assumptions on the mechanism of edema formation (detailed in Edema Formation section in S1 Text), our model showed local retinal thickening (Fig 11). A significant volume of fluid was observed 3 years post onset of initial capillary occlusion (Fig 11D). Extravascular fluid was responsible for the increase in retinal thickness. Moreover, the spot where fluid was present correlated with the boundary of the ischemic region. The current study is limited to qualitative illustration of how retinal thickening can be caused by abnormal VEGF synthesis by Mueller cells. This is the canonical function of VEGF producing leakage from capillaries. In the model edema is only a function of threshold VEGF levels and not a measure of local Starling type relationships. It is reflective within the model of loci of active leakage at sites of intact, non-occluded capillaries which also have elevated VEGF above a particular threshold. Edema within the model is a prediction of expected locations of either retinal thickening or of retinal cystic structures as seen clinically on ocular coherence tomography. The model does not address the actual range of cysts since these can often be present in areas of retinal ischemia lacking active leakage by fluorescein angiography [12] and thus can represent processes other than active leakage such as cellular necrosis or impaired retinal pigment epithelial pumping. The model does leave cysts present in areas in which capillaries had leaked but then subsequently become occluded similar to what occurs clinically.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. (CASE 1) retinal thickness over time.

Thickness of the retinal layer is represented by a color map at the end of year 0.5 (A), year 1 (B), year 2 (C) and year 3 (D). Color represents the magnitude in Z axis. The bluer the color, the thicker a local retinal area is. Fluid appeared between year 1 and year 2, at the venous edge of the area of occlusion. The flow network is overlaid upon the color map to present patent flow paths at the time point of observation. Color bar only represents the thickness of tissue but not the flow velocities. “FAZ” in the figure refers to foveal avascular zone, “A” in red refers to arteriole, and “V” refers to venule.

https://doi.org/10.1371/journal.pcbi.1004932.g011

Quantitative global measurement of disease progression

Following the onset of capillary closure, almost every additional occlusion caused average oxygen tension within cells to drop (Fig 12A). The fraction of hypoxic Mueller cells went up nearly linearly from week 72 to week 152 (Fig 12B). Distinct from the monotonic change in average oxygen tension and hypoxic fraction, the total volume inflow rate was non-monotonic, with continuous increases before a sharp decrease in week 124 and one more in week 152 (Fig 12C). Early stages of increased inflow rate might be attributed to physiological response to loss of a selected group of blood flow pathways. Because oxygen tension is assumed to be constant within inlet capillary blocks, as would be physiologically expected, a higher volume inflow rate gives more oxygen carried into the system within a given period of time. This might suggest that at early stages of the disease with few capillary closures, the system would compensate for the oxygen insufficiency by increasing the blood flow. In contrast, at later stages of disease when many capillaries connecting the arteriole to the venule were occluded, total flow declined and occlusions seemed to occur more frequently. To better map the model to clinically observed symptoms, our model used the scaled minimal cell-to-vessel distance d^min / a^MC as a metric to quantify disease progression spatially, where a^MC is the typical size of a Mueller cell (Fig 12D). This is basically the number of Mueller cell diameters to a patent, or unoccluded, vessel. The d^min / a^MC increased monotonically with each additional capillary closure. Interestingly, the capillary occlusion that triggered a rapid drop of the total volume inflow at week 124 had a mild effect on d^min / a^MC. This implied that topological location of a capillary segment was influential to the system blood flow supply, which might be neglected in a spatially-based metric. Thus, a flow-based model may be helpful in understanding the critical point at which the disease starts to evolve in a different way i.e. toward a decline in total inflow.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. (CASE 1) quantitative development of model properties with time.

(A) Average oxygen tension decreased with increasing occlusion events despite occasional small increases, and had an overwhelming drop near the end of simulation. (B) The hypoxic fraction of Mueller cells was observed to grow rapidly after the capillary occlusions induced by locally elevated VEGF, and near the end of simulation about one third of the Mueller cells were ischemic. (C) Total volume inflow rate established a rising trend before a major decline in week 124, after that it remained at a lower level than the starting point and eventually descended to below 90% of initial rate. (D) Average minimal cell-to-vessel distance maintained an increasing trend to reach 3.5 a^MC (Mueller cell diameter) units, which qualitatively reproduced the temporal pattern observed for hypoxic fraction of Mueller cells.

https://doi.org/10.1371/journal.pcbi.1004932.g012

The distribution of oxygen tension within all cells exhibited an essentially unimodal shape under the normal condition where most cells had oxygen tensions of 10 to 25 mmHg, a small portion of cells located near vessels had higher levels ranging from 35 to 40 mmHg (Fig 13) and no cells had an oxygen tension less than 4 mmHg O₂, i.e. no cells were ischemic. Capillary occlusions induced by elevation of VEGF gradually altered the distribution. An increasing number of cells turned hypoxic. The broad peak of cells at moderate levels of oxygen decreased and broadened with more cells both at lower oxygen levels and more cells from about 25–30 mmHg. This suggests that while the size of the hypoxic region in this tissue section was growing larger more cells were exposed to high oxygenation producing a bimodal oxygenation distribution.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 13. (CASE1) cellular oxygen distribution over time.

Fraction of cells in each 4 mmHg oxygen bin is shown for the normal condition and 3 additional times. The distribution of oxygen tension within all cells exhibited an essentially unimodal shape under the normal condition (normal-red bars) where most cells had oxygen tensions of 10 to 25 mmHg, a small portion of cells located near vessels had higher levels ranging from 35 to 40 mmHg and no cells had an oxygen tension less than 4 mmHg O₂. Capillary occlusions induced by VEGF gradually altered the distribution (week 0-blue bars, week 72-yellow bars, week 124-cyan bars). An increasing number of cells turned hypoxic. The broad peak of cells at moderate levels of oxygen decreased and broadened with more cells both at lower oxygen levels with each successive interval and more cells from about 25–30 mmHg in each successive interval. The cell oxygenation distribution gradually morphs from a unimodal distribution to a bimodal oxygenation distribution. Inset figure shows a comparison between normal condition and week 124 using a line-connecting-dot presentation, with finer oxygen tension spacing between two consecutive data points. The inset figure also strikingly shows the transition from unimodal to bimodal distribution as well as peak decreasing and broadening pattern at moderate levels.

https://doi.org/10.1371/journal.pcbi.1004932.g013

This model’s results are important both in terms of the images produced which bear a striking resemblance to those seen clinically and also as summarized in graphs showing changes occurring over time in a single run of the model for a specific initial capillary closure (Figs 8–13). Are these images and graphs consistent with what is known clinically about diabetic retinopathy? The pattern of capillary loss near the FAZ with expansion over time is very similar to what is seen in diabetic patients. The pattern of mixed ischemia and edema in the perifoveal area is also that usually seen clinically. The curve in Fig 12C is the most interesting in that flow rises from baseline levels for a period of time by as much as 13 percent and then declines. Clearly in end stage diabetic retinopathy with the entire capillary network occluded, the flow will go down. It is less clear that earlier stages in loss of the capillary network will result in increased total flow. This occurs because the model has some dilation of capillaries as a consequence of the adaptation module which would increase flows and also partially due to changes in network structure. The literature on blood flow in diabetes has been inconsistent due to measurements on diabetics at different stages of disease with a number of different technologies imaging flow in different locations [119]. The data from the Retinal Function Imager (RFI) measures flow velocities in small perifoveal vessels and seems most comparable to the vessel sizes in the model. This data is consistent with the model in that it shows an increased retinal blood flow velocity in diabetic patients without clinically seen morphological changes [119]. This would be similar to patients in the first year or two of the simulation (approximately 100 weeks).The percent increase in blood flow over controls was about 15% in [119], quantitatively similar to the model’s results. Physiologically this is possibly secondary to increased vasodilator mechanisms due to tissue hypoxia [120] and increased nitric oxide synthase [121] and possibly network changes which could have occurred even though the patients did not have clinically visible changes. The model replicates this effect though it does not explicitly utilize any analogous mechanisms. This pattern of increased macular blood flow accompanied by edema (Fig 11) is consistent with the regional distribution of diabetic lesions emphasized by [28]. Burgansky-Eliash 2010 has RFI data on patients with non-proliferative diabetic retinopathy [122] which show decreased blood flow at this stage of clinical retinopathy comparable to the blood flow decrease seen in the later weeks of this simulation. There is also clinical data on venous oxygen saturation in diabetes [123] showing that venous oxygen saturation in diabetes is elevated over that seen in normals. This model did not specifically treat venous oxygen levels but as blood flows were rising with a constant input oxygen saturation occurring simultaneously with a rise in hypoxic retinal cells, retinal oxygen extraction must be lower resulting in elevated venous oxygen saturation. The oxygen map shows this qualitatively as an expanded reddish tissue area (top right corner near “V” Fig 9C vs 9A) of elevated oxygen saturation around the venule. Hammer et al. (2009) interpret this elevation in venous oxygen levels to be due to a shortened arterio-venous passage time resulting in reduced oxygen extraction [123].

Replicate simulations

We performed 362 replications of the macular capillary sector simulation on Indiana University’s supercomputer BigRed II to pursue the consequences of different initial occlusion sites and explored the evolution of progression states in a flow-oxygen phase diagram (Fig 14). The probabilistic aspect of the capillary occlusions means that repeated runs of the model will not produce identical patterns of capillary loss but similarity of replications will be strongly influenced by network structure. In the flow-oxygen phase diagram, the normal condition (circles in dark blue) were clustered in a confined range, which showed that all simulations had similar equilibrium oxygen tension and total inflow rates initially as expected (Fig 14A). Circles gradually became scattered, because different simulations randomly picked different occlusion sites resulting in different disease progression trajectories. At 156 weeks, most simulations landed not far from the initial oxygen-flow states, which corresponded to situations without derived occlusions. On the contrary, some simulations showed distant end-point states situated mainly in low oxygen territory, which represented exacerbating progressive capillary occlusions. Of all such simulations, temporal trajectories were visualized (Fig 14B) if the simulations showed less than 75% total inflow rate or oxygen tension in year three. Most such simulations followed a clockwise trajectory temporally, starting from equilibrium state, transiting via low-oxygen and high-flow zone and ending in the low-oxygen and low-flow territory of the phase diagram, which indicated a severely damaged end-point of the capillary network. The rest seemed to pursue a temporal pattern of evolution but still remained in a transition stage. These graphs show the results of typical runs of the model. This might illustrate a general picture of how retinopathy progresses in terms of system oxygenation and total blood supply. Note that clinically many diabetics do not develop clinically significant retinopathy whereas a minority has significant propagation of capillary closure. This model would indicate that this variability of progression could have a significant element of probability, ‘bad luck’, in terms of which capillary was initially an occluded, in addition to other aspects of diabetic control and genetic/epigenetic factor [124–126]. It is now possible to gather time sequence data on the macular capillary network structure of diabetics through the utilization of AOSLO and in particular to do so in the two eyes of single patients, which possess different network anatomies but presumably identical environments in terms of glucose, blood pressure, and levels of systemic leukocyte activation. These data will allow comparisons of patient capillary occlusions and edema over time with sequentially constructed model based probability occlusion maps. This ought to allow refinement of model parameters or rejection of the model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 14. (CASE 1) Flow-Oxygen phase diagram of 362 replicate simulations.

Flow-Oxygen phase diagram describes how total volume inflow rate and oxygen tension evolve with time in all replicate simulations. Color represents simulation stages where bluish circles stand for early stages and reddish ones correspond to later stages. (A) Many simulations show clustering of flow-oxygen states through all monitored time points. Those clusters, quite close to the equilibrium state under normal condition, remained overlapped post the first occlusion. On the contrary, a few simulations displayed a scattered pattern of flow-oxygen states near the end of the simulation with either reduced flow or oxygen tension or both. (B) Temporal trajectories of the scattered points shown in A whose states had either less than 75% initial inflow rate or oxygen tension in year 3, were considered progressive. Simulations that lead to the most severe propagation of capillary occlusions tend to end with both lower than normal flow rates and lower oxygen tension.

https://doi.org/10.1371/journal.pcbi.1004932.g014

We further summarized the vulnerability of the capillary network given a certain initial occlusion site (Fig 15). Among all capillary segments carrying relatively slow blood flow, initial occlusion of two capillaries near the FAZ seemed to be most influential in triggering derived occlusions, while occlusions of others had a less significant impact. Both cases showed a spatially relevant patency distribution, with capillaries closer to the initial occlusion site bearing a higher frequency of occlusion. Not uncommonly, capillaries near the terminal FAZ venule and arteriole were also candidates of occlusion. Closure of these capillaries was likely to propagate injury to neighboring foveal arteriolar-venular sectors. These runs show similar results and show that certain portions of this subject’s capillary network are vulnerable to occlusion whereas others seem to be more resilient. Those more resilient areas tended to possess more densely situated capillaries and therefore the area of ischemia produced by a capillary occlusion is small and consequentially so is the resultant increase in local VEGF production and the probability of propagation. Occlusion or survival of a capillary is highly influenced by the local capillary network structure. This phenomenon seen in these results may be the result of more dense capillary networks further from the fovea where the visual impediment consequent to blood vessel opacity is less significant perhaps allowing more closely spaced capillaries. This pattern of loss of perifoveal capillaries, clinically called enlargement of the FAZ, is a commonly seen clinical pattern of macular ischemia development [11].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 15. Patency map of the macular capillary network.

The patency map illustrates the probability of each capillary segment remaining open given a specific initial capillary closure. These data are 2 subsets of the 362 replicate simulations in which 2 specific capillaries were occluded. (A) Patency map following CASE1 initial occlusion site shows that a high probability of secondary occlusions is limited to adjacent capillaries, largely those closer to the FAZ (n = 45). (B) Patency map following CASE2 initial occlusion site shows greater probability of broader propagation of injury (n = 30). Color represents the frequency of a capillary segment being patent after 3 years of simulated time since the initial occlusion. Warmer color corresponds to less vulnerability to occlusion, or higher capillary patency.

https://doi.org/10.1371/journal.pcbi.1004932.g015

Model improvements

Experimental data of a number of types can improve and validate the model. Arteriolar and capillary blood flow velocity data is needed to calibrate the network flow module and adjust the structural adaptation module. There are currently AOSLO based [127, 128] and ocular coherence tomography (OCT) angiography techniques [129] that will provide valuable inputs to future diabetic retinopathy modelling but no published techniques for measuring capillary blood velocities are sufficiently mature for routine use. Improved vascular imaging with capillary structures determined over the several retinal capillary layers in the macula would help extend the model beyond its current limitations. With regard to retinal edema, thickness projections can be validated and refined against OCT macular thickness maps and hydrostatic as well as oncotic pressures can be incorporated into the model. Refinements of the relationship of local retinal edema and capillary loss with visual functioning are also on the horizon [130] and can be assessed over time. More detailed structural modeling can be done of the cystic structure of the usual diabetic macular edema as can be determined by OCT. In the periphery the model predicts a gradual enlargement of the sorts of capillary free dark areas commonly seen on fluorescein angiography with occasional development of new small areas of capillary loss that also gradually expand and are blocked, at least temporarily by the oxygenated areas surrounding larger retinal vessels. Sequential fluorescein angiography, especially wide field angiography, would be validating in this area and allow refinement of the occlusion probability functions used in the model.

Comparison to a non-diabetic retinal vascular occlusive disease

Why does diabetic retinopathy progress to extensive areas of ischemia whereas certain other ischemic retinal vascular conditions do not progress? An example would be the arteriolar occlusions seen from particulates in IV drug users. The distinction between diabetic retinopathy and arteriolar occlusions caused by particulates is in the nature of the systemic inflammatory state in diabetes. This is likely the explanation for the non-progressive areas of capillary loss seen in certain retinal vascular conditions such as talc retinopathy in the IV drug user [138] in distinction to the enlarging areas of capillary occlusion in diabetic retinopathy. Diabetics have both a systemic inflammatory condition causing induction of cell adhesion molecules on leukocytes as well as induction of cell adhesion molecules throughout the retinal vasculature. The arteriolar occlusions seen in the IV drug user are local and likely do, through the creation of local ischemia, elevated local VEGF and induce ICAM receptors on local endothelial cells but do so without alterations in leukocyte cell receptors as the systemic inflammatory state is lacking. Therefore there is only local non-propagating retinal ischemia in the IV drug user and in animal embolic arteriolar occlusion models of diabetic retinopathy [139]. The absence of progressive capillary occlusion in either talc retinopathy or embolic models supports a ‘two-hit’ situation in diabetes, with induction of appropriate receptors on both the retinal endothelial cells and on the circulating leukocytes. Similar reasoning can be applied to the sharply geographically constrained ischemic damage produced by branch retinal vein occlusion.

The degree to which generalizations can be made to other organ systems is unclear but the presence of a positive feedback between the loss of a capillary and the likelihood of adjacent capillary loss may be less applicable to those organs which lack the constraint of a minimally sufficient, i.e. critical, vascular supply, imposed on the eye by evolution. This may partially explain the relative vulnerability of the eye in diabetes to the propagation of capillary occlusions and to the development of diabetic retinopathy.

A model provides a framework within which to analyze and interpret data. Currently, there are no appropriate models for either the progression of peripheral ischemia in diabetics or for diabetic macular disease. This model provides a structure within which to analyze and compare diabetic retinal vascular disease and some framework even if it proves incorrect or inadequate is required to move the understanding of this area forward in an era of refined imaging techniques able to provide inputs at the level required for capillary network modelling. This model, may, in some future instantiation, after considerable refinement and validation, allow improved treatment, with personalized prognosis, from imaging of a diabetic patient’s capillary network and with subsequent modelling allow a prediction regarding risk for vision impairment. Those at especially elevated risk may be candidates for earlier intervention either from anti-VEGF agents or from an alternative approach to photocoagulation. It seems reasonable that this sort of model could influence the course of patient therapy possibly in the not distant future.

Limitations

The current model is clearly limited relative to actual retinal vascular anatomy and physiology. The model was based on the advection and diffusion of the simplest known limiting factor, oxygen, rather than dealing with, say, removal of lactate or carbon dioxide. The model treats only one cytokine, a single VEGF isoform without consideration of other larger VEGF isoforms, PEDF, PDGF, erythropoietin, angiopoetin-1, angiopoetin-2 or angiopoetin-like 4, all of which likely play at least some role in vascular changes in diabetic retinopathy [117, 140–142]. In this case the level of therapeutic efficacy of the anti-VEGF injectable agents such as ranibizumab or aflibercept supports VEGF being a major [12, 100,143] if not the only important factor, supporting its relevancy and use in the model. This is in addition to all the evidence cited above (see Support for model’s physiological assumptions) supporting in choice of VEGF as the factor released from Mueller cells and ultimately responsible for capillary occlusion through elevation of ICAM and resultant leukostasis mediated capillary occlusion.

Movement of VEGF is a concern as there are a number of VEGF isoforms which have different molecular weights and also different structural domains which bind to tissue matrix [143]. Distribution of VEGF may depend on different states of the vitreous gel: attached, detached, syneretic, or absent (post-vitrectomy), as analogously, oxygen levels are dependent on intraocular location prior to, but not post, vitrectomy which presumably allows free advection [144]. We model only intra-retinal levels of VEGF, not vitreal levels, and there is little evidence specifically on these intraretinal levels and likely little relationship between these and the intra-vitreal levels seen late in the disease. It would seem that the current model is most applicable to the situation of non-proliferative disease with an attached intact vitreous resulting in minimal VEGF advection via fluid flow within the posterior chamber of the eye. This sort of elevated advecting VEGF would introduce another source of VEGF influence on local retinal areas not treated in the model.

The details of the complex mechanism of permanent capillary occlusion are unknown but it is likely that multiple intercellular adhesion events occur locally with brief capillary occlusions but also with endothelial cell damage and loss and that over time a local as well as a replacement endothelial cell population from the bone marrow are exhausted [145]. This seems to be another aspect of the impaired physiology of the diabetic. The model assumes a one-time irreversible occlusion event which we interpret as the local coup de grace for that capillary segment following a number of prior endothelial damaging leukocyte adhesion events. In the model occlusion probability was assigned a formula related to local VEGF levels as well as vessel size and flow. This is a sensible approximation to a complicated recurrent local adhesion/occlusion event but has no specific supporting data from the literature.

There are also limitations in the model’s approximation to the capillary anatomy in the macula. The model is three dimensional but the capillary bed treated is only two dimensional in the sense that it is a single lamina thought of as the retinal ganglion cell layer capillary network. It is a single capillary layer but only for about 60 microns from the FAZ. There are currently technological difficulties with resolution of the several capillary layers present further from the FAZ and in this area the model is based only on the inner capillary network imaged by AOSLO. The macular thickness increases over the clivus and thus the amount of tissue dependent on the capillary bed is changing with distance from the FAZ.

A number of clinical signs of diabetic retinopathy were not modelled such as microaneurysms, intra-retinal microvascular abnormalities, cotton wool spots, or dot/blot hemorrhages. These signs can reflect local ischemia and locally elevated VEGF but we feel they are not etiologically related to the topic of progression of capillary occlusions.

The word ‘periphery’ is used loosely in ophthalmology. A model limitation regarding the retinal periphery is that distinct and varying retinal vascular anatomy is seen from the most posterior retinal periphery, that near the arcades (Fig 1), to that anterior to the eye’s equator and especially that most anterior retina near the ora serrata (S8 Fig). We have modelled the extreme periphery (anatomically a single capillary layer) and made statements about ‘peripheral retina’ in general, much of which posterior to the equator has 2 capillary layers.

Summary

We have implemented a quantitative model of the impact of diabetes on progressive capillary occlusion in retinal capillary networks producing a number of qualitative results very similar to the clinical picture in diabetic retinopathy. A full vascular model from a retinal arteriole to a venule with the linking capillaries is included with oxygen advection, oxygen diffusion, and oxygen consumption. This model requires only the capillary network anatomy, arteriolar and venular pressures, and a number of physical parameters such as a diffusion constant for VEGF and for oxygen. The model produces capillary occlusion on a stochastic basis dependent on local VEGF levels and flow where the local VEGF level is seen physiologically as elevating ICAMs, increasing leukocyte adhesion and ultimately resulting in the capillary occlusion. This results in a local adverse feedback cycle often producing progressive capillary occlusions and enlarging contiguous areas of retinal ischemia. With these simple properties if the model is implemented using a capillary network for a perifoveal area it shows regions of ischemia and macular edema similar to what is seen clinically in the macula including enlargement of the foveal avascular zone. If implemented on a vascular pattern similar to that seen in peripheral retina, model damage proceeds more in line with what is seen clinically in the periphery developing large ischemic areas without edema. This model seems to show that the type and nature of retinal damage depends largely on the pre-existing capillary network morphology, suggesting that it may, after further refinement and validation, have value not only by providing a model through which to augment our understanding of diabetic retinopathy but even potentially, after sufficient clinical validation, for prognostic evaluations and timing of treatment in individual patients.