Towards a Glaucoma Risk Index Based on Simulated Hemodynamics from Fundus Images

Abstract

Glaucoma is the leading cause of irreversible but preventable blindness in the world. Its major treatable risk factor is the intra-ocular pressure, although other biomarkers are being explored to improve the understanding of the pathophysiology of the disease. It has been recently observed that glaucoma induces changes in the ocular hemodynamics. However, its effects on the functional behavior of the retinal arterioles have not been studied yet. In this paper we propose a first approach for characterizing those changes using computational hemodynamics. The retinal blood flow is simulated using a 0D model for a steady, incompressible non Newtonian fluid in rigid domains. The simulation is performed on patient-specific arterial trees extracted from fundus images. We also propose a novel feature representation technique to comprise the outcomes of the simulation stage into a fixed length feature vector that can be used for classification studies. Our experiments on a new database of fundus images show that our approach is able to capture representative changes in the hemodynamics of glaucomatous patients. Code and data are publicly available in https://ignaciorlando.github.io.

The original version of this chapter was revised: An acknowledgement has been updated. The correction to this chapter is available at https://doi.org/10.1007/978-3-030-00934-2_106

1 Introduction

Glaucoma is a neurodegenerative condition that is the leading cause of irreversible but preventable blindness [1]. The disease is characterized by the interruption in the communication between the retinal photoreceptors and the brain, preventing visual signals to be carried to the brain via the optic nerve [1]. These progressive changes are asymptomatic until advanced stages of the disease, when undiagnosed patients notice the disease when a considerable amount of visual field has been lost. The major treatable glaucoma risk factor is the intra-ocular pressure (IOP), although individuals with high IOP might not develop glaucoma, and subjects with normal IOP values can develop it [2]. Hence, alternative genetic, vascular, anatomical or systemic factors are being extensively studied to better understand the pathophysiology of glaucoma [2, 3].

The ocular hemodynamics has long been observed to be affected by glaucoma [2, 4] through imaging modalities such as color Doppler imaging or laser Doppler flowmetry, among others [4]. These tools provide non-invasive measurements of hemodynamic parameters from ocular vascular structures such as the central retinal artery (CRA), the ophthalmic artery (OA) or the retinal capillary network. In turn, the analysis of the hemodynamic characteristics of the retinal arterioles and venules, easily accessible through fundus photography, has been limited despite the availability of this imaging modality for several decades.

In this paper, we propose a first approach to characterize the hemodynamics of the retinal microvasculature of glaucomatous and control subjects based on computer simulations performed on vascular networks extracted from fundus images (Fig. 1). Our method is based on patient-specific graph representations of the arterial network that are used as the topological vascular substrate to build 0D models, which account for the steady, incompressible flow of a non Newtonian fluid (the blood) in rigid domains. Although simpler than fully 3D [5] or 2D [6] blood flow models, 0D models (also known as lumped parameter models) are computationally cheaper and, thus, enable the simulation of blood flow in large networks of vessels [7]. Similar modeling approaches have been used before in the context of retinal hemodynamics simulation. In [7], the blood hematocrit is incorporated as transport species in the model, considering constant cross-sectional radius per vessel segment in a network of arteries and veins of a mouse model. More recently, 0D modeling has been applied in [8] for characterizing abnormal changes due to diabetic retinopathy. Although the simulation approach is similar to ours, a Newtonian fluid model was used in [8]. Furthermore, instead of using the whole retinal network, the numerical experiments have been performed in a single vascular bifurcation, which was manually isolated for each image. As a further novel aspect of the present study, the simulation outcomes are summarized in a fixed-length feature vector using a novel technique inspired in the bag of words [9] technique for feature representation. This allows the characterization approach to be applied in combination with linear classifiers to study its relationship with the presence/absence of the disease. We empirically validate this approach on a new database, namely LES-AV, comprising 22 fundus photographs with available manual segmentations of the retinal vessels and their expert classification into arteries and veins.

In summary, the key contributions of this study are: (i) a computational workflow for quantifying patient-specific hemodynamics non-invasively; (ii) a feature representation technique to comprise the simulated outcomes into a fixed length feature vector; (iii) a first pilot study showing that glaucomatous patients exhibit changes in the hemodynamic variables with respect to control subjects; and (iv) a new database of fundus images that is released, jointly with code and implementation details, in https://ignaciorlando.github.io.

Fig. 1.

Proposed method for characterizing patient-specific retinal hemodynamics.

2 Materials and Methods

The experiments were conducted on a new database of color fundus images, called LES-AV, comprising 22 images/patients with resolutions of $30^\circ$ FOV and $1444 \times 1620$ pixels (21 images), and $45^\circ$ FOV and $1958 \times 2196$ pixels (one image), with each pixel $=$ 6 $\upmu$m. Demographic data about the groups are provided in Sect. 3 and the supplementary materials.

2.1 Preprocessing

Patient-specific segmentations of the retinal vasculature must be first retrieved from a fundus image, and subsequently classified into arteries and veins. Although several methods have been introduced for this purpose [10], we follow a semiautomatic approach to ensure a proper input for the simulation step. Hence, a first coarse segmentation of the vessels is obtained using a fully convolutional neural network [11] (see supplementary materials for details). The segmentation is then manually refined to improve the vessel profiles and their connectivity, and each structure is manually labelled as artery or vein. Subsequently, the arteries are taken and their corresponding centerlines are automatically identified using a skeletonization method [12]. The result is a binary mask T with $N_I$ connected components. Since the vasculature is observed as a 2D projection of the original 3D structure in a fundus image, the vessels in the optic disc (OD) are generally overlapped with each other, making unfeasible the identification of the roots of the trees in any realistic scenario. To mitigate this issue, the arterial segments inside the OD are removed from T, and a root pixel $I_t$ is automatically chosen for each connected component $T_t \in T$, such that $I_t$ corresponds to the closest pixel to the OD center. Subsequently, each tree $T_t$ is modeled as a graph structure $G_t = \langle I_t, \mathcal {V}_t, \mathcal {E}_t \rangle$, where each node $v_j \in \mathcal {V}_t$ is a set of branching pixels in $T_t$, and each edge $e_{(i,j)} \in \mathcal {E}_t$ is a segment of centerline pixels connecting two branching points $v_i$ and $v_j$. The pixels in the skeletonized arterial trees are automatically classified as part of a node or an edge in $G_t$ by analyzing their neighborhood. Finally, the cross-sectional lumen radius at each centerline pixel is estimated from the segmentation using the Euclidean transform [13].

2.2 Simulation of the Retinal Hemodynamics

The simulation is performed on a patient-specific computational mesh associated to the graphs $G_t$ (Sect. 2.1). This structure comprises the key information required by the model to perform the simulation: (a) the topological connectivity of the arterial segments, through their branching points; (b) the pixel-wise information of the arterial radius, $r_i$, $i=\{1,...,N\}$, with N the number of centerline pixels (a computational node per pixel is considered); (c) the root points, which are the network inlets $I_t$, $t=1,...,N_I$; and the network outlets $O_m$, $m=1,...,N_O$.

A steady state 0D model is used to describe the incompressible flow of a non-Newtonian fluid in non-compliant vessels. Such a model is ideal for this scenario, as it allows to efficiently analyze different simulations, and provides physiologically reasonable coarse descriptions of the hemodynamics in vascular networks. The non Newtonian properties of the blood are incorporated by expressing the viscosity as a function of the vessel cross-sectional radius, to describe the Fahraeus-Lindqvist effect [14]. As a result of the simulation process, the blood pressure ($P_i$) and the blood flow rate $Q_i$ are computed at each computational node (i). The governing 0D equations are detailed below.

The standard mass conservation and pressure continuity junction model is used at each bifurcation point, where node i branches to nodes j and k (Eq. (1)).

Given an arterial segment of M pixels length, $M-1$ lumped parameter elements are used, where the mass conservation and the hydraulic analogue of the Ohm’s law are considered. Then, for an element formed by nodes i and $i+1$ (Eq. (2)), the resistance to the flow at each segment is modeled with the Poiseuille’s law: $R_{i,i+1} = 8\mu L / \pi r_{i,i+1}^4$ where L is the segment length, $r_{i,i+1}=(r_i+r_{i+1})/2$ is the average radius and $\mu =\mu (r_{i,i+1})$ is the blood viscosity, which is assumed to be a function of the arterial radius (since $r<150$ $\upmu$m), see the supplementary materials.

$$\begin{aligned} {\left\{ \begin{array}{ll} Q_i &{}= Q_j + Q_k,\\ P_i &{}= P_j = P_k. \end{array}\right. } \end{aligned}$$

(1)

$$\begin{aligned} {\left\{ \begin{array}{ll} Q_{i+1} &{}= Q_i, \\ Q_{i+1} &{}= \dfrac{P_i - P_{i+1}}{R_{i,i+1}}. \end{array}\right. } \end{aligned}$$

(2)

Regarding boundary conditions, the pressure and flow rates, are prescribed at the inlet and outlets, respectively:

$$\begin{aligned} {\left\{ \begin{array}{ll} P_{I_l} = P_0, &{} \forall l=1,\ldots ,N_I, \\ Q_{O_m} = \beta r_{O_m}^\gamma , &{} \forall m=1,\ldots ,N_O. \end{array}\right. } \end{aligned}$$

(3)

The pressure at all the inlets is set to the mean CRA pressure, $P_0$, while the flow at the outlets is set to the value given by Murray’s law, which relates the flow rate to the outlet radius to the power of $\gamma$. Specifically, we set $\gamma =2.66$ [15], and putting $\beta = Q_T / \sum \limits _m r^\gamma _{O_m}$ ensures the total retinal flow $Q_T$.

The inlet pressure is set to the mean value of normotensive patients (i.e. $P_0=62.22$ mmHg). Due to the lack of consensus on normal values for $Q_T$  [16], we proposed three scenarios, namely $Q_T^{\text {SC1}}=30, Q_T^{\text {SC2}}=45.6, Q_T^{\text {SC3}}=80\,\upmu$l/m, which represent the mean (SC2) and extreme values (SC1 and SC3) for healthy subjects. These values are preset for any given subject, without using patient-specific parameters. Further details about the physiological considerations of the parameters choice are provided in the supplementary materials.

The resulting real linear system of equations is solved using LAPACK dgles function through a QR factorization within a custom C++ source code. After solving the system, for each centerline pixel i in the network G, a feature vector $\mathbf {F}_i$ is computed, comprising: the flow rate $Q_i$, the blood pressure $P_i$, the mean cross sectional velocity $v_i=Q_i/(\pi r_i^2)$, the Poiseuille’s resistance $R_i$, the Reynolds number $\text {Re}_i=\rho 2 r_i v_i/\mu _i$ and the wall shear stress WSS$_i=4\mu _i Q_i / (\pi r_i^3)$. Here, $\rho =1.040$ g/cm$^3$ is the blood density. At bifurcation points, the variables assume the value of the last computational node of the parent artery.

2.3 Bag of Hemodynamic Features (BoHF)

One key aspect of the functional characterization is how to construct a feature vector $\mathbf {F}$ for discriminating control subjects from diseased. As $\mathbf {F}_i$ is given for every centerline pixel , it is necessary to find an alternative feature representation with a fixed length, suitable e.g. for linear discriminant analysis. To this end, we propose a new technique inspired by the Bag of Words (BoW) method [9], named as Bag of Hemodynamic Features (BoHF).

Our analysis is focused on three relevant structures for the graph $G^{(u)}$ of the u-th subject: the terminal and the bifurcation points, and the arterial segments. For every pixel i being a terminal or a bifurcation point, the vector $\mathbf {F}_i$ is used. However, for a given arterial segment $s \in \mathcal {E}$, the average vector $\tilde{\mathbf {F}}_s$ over the pixels enclosed by the segment is used as a summary statistic. This allows our representation to lower its variance and to normalize for the length of an arterial segment. The resulting vectors associated to each arterial segment, bifurcation and terminal points are collected into a set, denoted X.

Let $S = \{ (X^{(u)}, y^{(u)}) \}, u \in \{1, ..., n\}$ be a training set of n samples, where $X^{(u)}$ is the set of summary statistics at arterial segments, bifurcation and terminal points described above, and $y^{(u)} \in \mathcal {L} = \{ -1, +1 \}$ is its associated binary label (e.g. healthy or glaucomatous). Since the size of $X^{(u)}$ varies from one subject to another, our purpose is to map it into a feature vector $\mathbf {x}^{(u)} \in \mathbb {R}^d$, with d fixed. In the traditional BoW approach, the $\mathbf {x}$ vectors are histograms of codewords (analogous to words in text documents), and are obtained based on a codebook $\mathcal {C}$ (analogous to a word dictionary). $\mathcal {C}$ must be designed such that it summarizes the most representative characteristics of each class in $\mathcal {L}$. Let us denote by $S_{-1}$ and $S_{+1}$ the set of negative and positive samples in S, respectively. A codebook $\mathcal {C} = \{ \mathcal {C}_{-1}, \mathcal {C}_{+1} \}$ can be automatically learned by applying k-means clustering in the p-dimensional feature space of each subset $S_{(\cdot )}$. The centroids of each of the learned k clusters can be then taken as the codewords, resulting in a set of k codewords for each class ($|\mathcal {C}| = d = 2k$). Hence, the feature vector $\mathbf {x}$ can be obtained such that the j-th position in $\mathbf {x}$ corresponds to the number of nearest neighbors in $X^{(u)}$ to the j-th code in $\mathcal {C}$. If the codewords are representative of each label group, it is expected that the samples belonging to the positive/negative class will have higher/lower values in the last k positions of $\mathbf {x}$.

Finally, the discrimination of glaucomatous patients from control subjects can be posed as a binary classification problem that can be tackled using any binary classification algorithm, e.g. $\ell _2$ regularized logistic regression. The linear discriminant function of the logistic regression model $f(\mathbf {x}) = \langle \beta , \mathbf {x} \rangle$ with parameters $\beta$ is obtained by solving: $\hat{\beta } = \mathop {\mathrm {arg\,min}}\nolimits _\beta \bigg [\lambda \Vert \beta \Vert _2^2 + \sum _{u=0}^{n} \log (1 + e^{-y^{(u)} \langle \beta , \mathbf {x}^{(u)} \rangle })\bigg ]$.

Fig. 2.

Mean radius vs. flow rate per segment for healthy (black) and glaucomatous (red) subjects in the LES-AV database. (a) Samples per image and fitted curves for each individual. (b) Fitted curves for each group. The SC1 and SC3 were used for the glaucomatous and healthy groups, respectively.

3 Results

The computational simulations of all scenarios resulted in mean values of pressure drop ($\varDelta P$ to inlet, $P_0$) and blood velocity (v) consistent with previous studies in the human retinal circulation [6].

The ability of the BoHF to capture the hemodynamic changes in the glaucomatous (G) from normal (H) subjects was empirically assessed using $\ell _2$ regularized logistic regression on LES-AV, in a leave-one-out cross-validation setting. The SC2 simulations were used for all the cases to avoid any bias in the classification performance. The values of $k \in \{2, ..., 15\}$ (BoHF length) and $\lambda \in 10^i, i \in \{-3, ..., 0, ..., 5\}$ (regularization parameter) were fixed to the values that maximized the accuracy in a held-out validation set randomly sampled from the training set. The features were centered to zero mean and unit variance using the training sets mean and standard deviation. The resulting area under the ROC curve was 0.70248, with an accuracy of 68.18%.

The hemodynamic quantities obtained using the outcomes of SC1 and SC3 for the G and H groups, respectively, are in line with clinical observations. Figure 2(a) shows the flow vs. mean radius per segment for all subjects. The exponential curves fitted for each group are depicted in Fig. 2(b). The variables are correlated (Spearman’s $\rho =0.71$, $p\ll 0.01$), which is consistent with the literature [6].

Finally, statistical analysis at a confidence level of 99% were performed on a per-patient and a per-measurement basis. On the former (Table 1), we focused on the age, sex and number of segments (NoS). The differences in the mean values of the NoS and the age were not significant (Wilcoxon test, $p>0.1$). Moreover, the patient diagnosis was independent from the sex ($\chi ^2$ test, $p=1$).

Table 1. Summary of statistics at a per-patient level. None of the variables present statistically significant differences between groups ($p\gg 0.01$).

Table 2. Summary of statistics at a per-measurement level (terminals, bifurcations and arterial segments). ($^*$) points statistically significant differences between groups ($p\ll 0.01$). ($^\dagger$) reached $p=0.02$.

On a per-measurement basis (Table 2), all the arterial segments, terminal and bifurcations points were considered as separate measurements. The analysis was focused on the $\varDelta P$, v, r, age and sex. The differences in the mean value of $\varDelta P$ and v were larger in the H group, while the age was larger in the G group ($p\ll 0.01$) (two tailed Wilcoxon test). The r was lower in the G group, although not significant ($p=0.31$). As in the per-patient analysis, the diagnostic was independent from the sex ($\chi ^2$ test, $p=0.02$). The Pearson correlation coefficient between age vs. $\varDelta P$ and age vs. v was close to zero and not significant, suggesting that the age is not an influential factor for the comparison of the means.

4 Discussion

We have presented a first approach to characterize the retinal hemodynamics of glaucoma patients using computational fluid-dynamics. To the best of our knowledge, this is the first study in which the relationship between glaucoma and simulated hemodynamic outcomes of the retinal arterioles are studied.

The arterial radius was previously observed to be smaller in glaucoma patients than in normal subjects [17], a setting that is consistent with our data. On the other hand, in-vivo velocity measurements in the CRA of affected patients was shown to be lower than in control subjects [18]. Under these considerations, we have used SC1 for the G group and SC3 for the H group in our first pilot study, resulting in a simulated behavior that is in line with existing observations made on other measurable vessels (i.e. $\text {G}(v)<\text {H}(v)$) [18]. Moreover, this rendered physiologically correct values of $\varDelta p$ in the retinal network.

The classification experiment based on BoHF and the unbiased SC2 outcomes showed an AUC different than random for identifying the glaucomatous patients, indicating that the tool is able to preserve the variations in the hemodynamics of the arteries affected by the disease. Future efforts will be focused on incorporating the venous network and also modeling of the IOP, which is relevant for clinical studies [19]. Also, the incorporation of a deep learning based method for simultaneous segmentation and classification of the vasculature will be explored in order to allow a more efficient characterization of larger populations.

Finally, our MATLAB/C++/python code and the LES-AV database are publicly released. To the best of our knowledge, our data set is the first in providing not only the segmentations of the arterio-venous structures but also diagnostics and clinical parameters at an image level.