In-silico model of hypoxic areas caused by micro-embolisms in the rat brain

(1)

In-silico model of hypoxic areas caused by

micro-embolisms in the rat brain

Sjoerd Terpstra

Report Bachelor Project Physics and Astronomy

Conducted between February 8, 2021 & May 31, 2021

Student number 11251980

Universities University of Amsterdam (UvA) & Vrije Universiteit (VU) Study programme Physics & Astronomy (BSc JD)

Course Bachelor Project Physics and Astronomy (15 EC) Institute AMC Biomedical Engineering and Physics

Examiner Prof. A.G.J.M. van Leeuwen Supervisor Prof. E.T. van Bavel

Daily supervisors MSc T. Georgakopoulou & MSc N. Arrarte Terreros Submission date June 1, 2021

(2)

Abstract

Endovascular treatment is commonly used to treat patients with acute ischaemic stroke. However, around 50% of the patients still remain functionally dependent after successful recanalization. It is hypothesised that this is caused by micro-thrombi, released from the original thrombus during endovascular treatment, which are transported to the distal arterioles in the brain. These can cause numerous small hypoxic regions. Georgakopoulou et al. (2021) showed that microspheres, similar to micro-thrombi, can create numerous hypoxic areas in rat brains. In the current study, a stochastic model that generates hypoxic areas caused by microspheres according to a Gaussian risk function is proposed. This study aims to reproduce the spatial relationship between microspheres and hypoxic regions using data from Georgakopoulou et al., 2021. In addition, this study investigates the cumulative effect of microspheres in causing hypoxic areas. The model approximately reproduces the spatial relationships between microspheres and hypoxic areas and shows that at least two microspheres need to cooperate to create hypoxic areas.

Samenvatting

Een acute ischemische beroerte wordt veroorzaakt door een bloedprop in een van de grote bloedvaten van de hersenen. Om deze patiënten te helpen kan de bloedprop verwijderd worden door met behulp van een stent retriever via de slagaderen de bloedprop te pakken en uit het lichaam te halen. Hierdoor kan het bloed weer gaan stromen en de patiënt herstellen. Helaas blijft rond de helft van de patiënten afhankelijk van de zorg van anderen na deze ingreep. Een hypothese om dit verschil tussen een succesvolle ingreep en een slechte uitkomst voor de patiënt te verklaren, is dat door deze ingreep vele kleine bloedpropjes ter grootte van een paar micrometer losraken van de grote bloedprop. Deze propjes worden naar de haarvaten in de hersenen getransporteerd, waar ze de bloedtoevoer kunnen blokkeren. Hierdoor heeft het weefsel in de omgeving van de blokkades zuurstof tekort en kunnen de cellen in het weefsel afsterven. Dit wordt ook wel hypoxisch weefsel genoemd. Omdat de propjes zo klein zijn, is het met de huidige beeldvormende technieken niet mogelijk dit proces direct in patiënten te bestuderen. Daarom zijn er experimenten gedaan op ratten om de relatie tussen deze propjes en hypoxisch weefsel in kaart te brengen. De data van deze experimenten laten zien dat de propjes inderdaad hypoxisch weefsel kunnen veroorzaken, en dat er een ruimtelijke relatie is tussen de twee. Om dit beter in kaart te brengen hebben we een computer model gemaakt dat aan de hand van de positie en grootte van de propjes, de positie en grootte van het hypoxisch weefsel tot op een redelijke hoogte kan voorspellen. Het blijkt dat meerdere propjes nodig zijn om samen hypoxisch weefsel te veroorzaken.

(3)

Acknowledgements

I want to thank my supervisor Ed van Bavel for his time to help me guide this project. The meetings about my project were always very insightful and I value the discussions we had on this topic. I also want to thank my daily supervisors Sissy Georgakopoulou and Nerea Arrarte Terreros for their input. Thank you for your time to help me find my way in the previous research and data for this topic, and for giving me feedback on the thesis. Furthermore, I would like to thank my examiner Ton van Leeuwen. Thank you for your interest in my thesis. Finally, I want to thank the research group at the Biomedical Engineering and Physics department. It was really interesting to witness the research done in the department.

(4)

1 Introduction

Stroke is one of the leading causes of death worldwide, ranking only second behind ischaemic heart disease (World Health Organization, 2020). Acute Ischaemic Stroke (AIS) is the most prevalent type of stroke, representing up to 80% of all cases, and unfortunately up to 50% of survivors of AIS remain functionally dependent on others after a stroke (Donkor, 2018).

According to Phipps and Cronin (2020, p.1), “(...) [AIS] is defined by the sudden loss of blood flow to an area of the brain with the resulting loss of neurological function”. AIS is caused by the occlusion of a cerebral vessel due to an embolism, a blockage in a blood vessel caused by a thrombus. This deprives brain areas of the needed blood flow and consequently reduces the oxygen supply to the neural tissue. This can lead to irreversible damage and thus loss of neurological function if the vessel remains occluded over a long period of time. That is why it is so crucial to treat patients in the first couple of hours after the onset of symptoms of AIS. There are two main treatment options for patients: thrombolysis and endovascular treatment (EVT), which can either be used individually or combined.

Thrombolysis is a treatment where the thrombus is dissolved using drugs. The most common drug used is recombinant tissue plasminogen activator (rt-PA). Although this drug can improve the recovery of patients, it has serious side effects. rt-PA can cause severe bleeding in the brain, which can often lead to death (Wardlaw et al., 2014). Therefore, efforts are made to develop alternative treatment options that preferably have less serious side effects and/or improved clinical outcome. One such alternative is EVT.

EVT is a relatively new treatment where the thrombus is mechanically removed using mechanical thrombectomy (Bhaskar et al., 2018). During this procedure, a stent retrieving device is inserted into the cerebral vessel at the site of the occlusion, from an accessible blood vessel. The stent traps the thrombus, as Figure 1 shows. The device then retracts it out of the occluded vessel, hopefully resulting in successful recanalization of the vessel.

EVT improves the health outcome of AIS patients over other treatment options (Goyal et al., 2015), but still a majority do not (fully) recover. One hypothesis that explains a part of the discrepancy between successful recanalization and the complete recovery of patients, is that the procedure can lead to small thrombus fragments being released from the thrombi into the cerebral blood flow. These micro-thrombi are then transported to the distal capillaries of the brain (Chueh et al., 2016). These can cause numerous micro-embolisms, leading to lower perfusion of oxygen in multiple parts of the brain, resulting in numerous small ischaemic areas, infarcts and hypoxic areas.

Figure 1: During endovascular treatment of acute ischaemic stroke, a stent is inserted in the occluded vessel. The stent envelops the thrombus and then retracts it out of the vessel. Source: INSIST, 2021.

(6)

1 INTRODUCTION

cells to survive over a longer period of time. There is no hard cut-off for the level of oxygen at which tissue gets hypoxic. Instead, this is defined by the level for which the cells in the tissue cannot function anymore, which can differ between regions in the body and between individuals (Samuel and Franklin, 2008). In this study HYP are defined as tissue with pO2 < 10 mmHg. Since HYP can lead to cell death, it can play a very important role in

the functional outcome of patients after EVT treatment.

To successfully combat AIS, the previously mentioned treatment options have to be continuously improved. Since clinical trials are both expensive and difficult to perform, because of technical and ethical restrictions, in-silico models are used to improve and judge treatments before clinical trials are performed. The aim of the INSIST (IN-Silico trials for treatment of acute Ischemic STroke) project1 _{is to build a network of models}

to test and improve treatment options for AIS (Konduri et al., 2020). This is achieved by creating virtual populations of patients, and by applying models to simulate AIS and virtual treatments to these patients to measure their outcomes. In order to create these models, various aspects of AIS need to be further examined to be able to accurately predict the outcome of patients based on treatments. The possible negative effect of micro-emboli is one of them.

Unfortunately, it is not feasible to track micro-emboli in-vivo with current imaging techniques, due to their limited resolution (Chueh et al., 2016; van Veluw et al., 2017). Although progress is being made to detect smaller lesions, most imaging techniques have resolutions on a scale of mm’s and not on µm’s (van Veluw et al., 2017). This makes it difficult to directly study the possible effect of micro-emboli on the health outcome of patients. To this end, in-silico models and rodent models are used instead.

Georgakopoulou et al. (2021) performed an in-vivo study on rats to examine the influence of micro-emboli on the perfusion and oxygenation in the brain. They injected polystyrene microspheres (MS) of sizes 15, 25 and 50 µm via the common carotid artery in the brains of 18 rats. After each of 1, 3, and 7 days, six animals were sacrificed and slices of 500 µm of their brains were imaged. In this study they found that MS can cause multiple ischaemic areas, infarcts and HYP.

In Georgakopoulou et al. (under review), the spatial relationships between MS and HYP from the study of Georgakopoulou et al. (2021) were quantified using G and K functions. These functions describe the clustering and spread of events in relation to each other (O’Sullivan and Unwin, 2010). This study found that both the MS and HYP were clustered. Also, a correlation between the MS and HYP were found up to a distance of 1000 to 1500 µm. The findings imply that multiple MS can work together to create HYP in the brain.

The exact mechanism that underlies this relationship still remains unclear. In this study, an in-silico model is proposed to investigate this process further. The model is used to simulate the HYP caused by MS. The data for the MS and HYP are taken from the rats that were sacrificed at day 1 from Georgakopoulou et al. (2021). The model computes a risk matrix of the brain, where each MS causes a risk area based on a Gaussian distribution. This risk matrix is segmented to a risk image, where the domain is divided between HYP and non-HYP on the basis of a threshold.

The first objective of this model, referred to as the Gaussian risk model (GRM), is to examine if the quantitative spatial relationships between MS and HYP from Geor-gakopoulou et al. (under review) can be reproduced in this model. The second objective is to further investigate how the MS work together to form HYP.

(7)

1 INTRODUCTION

The motivation of this study is to better understand how MS cause HYP. Hopefully, it can shed more light on the interactions between MS, and their influence on hypoxic regions. Focusing on the effects of micro-emboli caused by EVT is informative in the effort to improve clinical outcome of patients.

(8)

2 GAUSSIAN RISK MODEL 2.1 Data input

2 Gaussian risk model

The in-silico model that is developed is a stochastic model with as aim to simulate HYP based on data of MS. The general idea is that each MS creates an area in the brain at risk for becoming hypoxic. This risk area is based on a Gaussian distribution, hence the model will be referred to as the Gaussian risk model (GRM); a model that models the risk of HYP based on Gaussian distribution functions.

The GRM consists of 4 main parts:

1. Data input. The first part handles the input data for the model, as partially summarised in Table 1. In addition to this, the centroids of both the MS and HYP are supplied to the model.

2. Generating the risk matrix. In this part the risk matrix is calculated. This matrix denotes the risk of having a hypoxic region at each point in a discretized rectangular cuboid of the brain.

3. Risk matrix segmentation. In this part the risk matrix is segmented using a threshold to a risk image. This risk image contains the hypoxic regions.

4. Statistics of output. The final part contains information about the output of the simulation and the calculation of the metrics used to judge the output.

Each of these modules are explained in the following sections.

2.1 Data input

Table 1: The input data for the GRM per animal. The number of total MS and the number of MS per diameter is shown. The total number of HYP and the total HYP volume are also shown. GRM: Gaussian risk model, HYP: hypoxia, MS: microspheres.

Animal Total MS 15 µm MS 25 µm MS 50 µm MS Total HYP Volume HYP (µm3) 1 187 158 27 2 152 2.78 × 108 2 78 61 17 0 83 9.12 × 107 3 41 37 4 0 27 2.36 × 107 4 102 102 0 0 70 6.68 × 107 5 112 100 12 0 67 1.79 × 108 6 102 77 23 2 94 1.47 × 108

In this study, all data used are collected and processed by Georgakopoulou et al. (2021). Only the data of the animals who were sacrificed at day 1 are used in this study, because it was shown by Georgakopoulou et al., 2021 that some HYP recover after day 3 and day 7. This mechanism is not part of this study, since the focus here is on the causal relationship between MS and HYP, and not on the recovery mechanism.

The data of each animal consists of information about the MS and HYP in a brain slice of 500 µm. The number, diameters and positions of the MS, and the number, volumes and centroids of the HYP are processed. Table 1 shows a summary of these data.

To generate the risk matrix, the GRM uses the number, diameters and positions of the MS. The total volume of the HYP is used to segment the risk matrix to a risk image

(9)

2 GAUSSIAN RISK MODEL 2.2 Generating the risk matrix

Figure 2: An example of the MS input for the GRM. This is a projection of all MS on the (x, y)-plane. The MS are the white dots, artificially scaled to 40x40 µm to make them better visible for the eye. GRM: Gaussian risk model, MS: microspheres.

(see subsection 2.3), to obtain the HYP generated by the model. The data of the number, volumes and centroids of the HYP are compared to the model output.

Figure 2 shows the MS input for the first of the six animals as a projection on the (x, y)-plane of a rectangular cuboid domain (a 3D equivalent of a rectangle). The white dots denote MS, all scaled to 40x40 µm for this plot to ensure good visibility. It is worth noting that the rectangular cuboid domain encompasses the whole segmented part of the brain and thus all MS, but since the brain has a different shape than a rectangle, the domain is larger than the brain slice. Therefore, in the corners of the domain, where there is no brain, there are no MS found. In this study however, there is no distinction made between areas that are, and are not part of the actual brain. This makes it possible for HYP to develop in parts outside the brain boundaries, but inside the rectangular cuboid domain.

2.2 Generating the risk matrix

The brain volume is discretized by creating a 3D risk matrix Mrisk, whereby each element

represents a voxel in a rectangular cuboid space with coordinates (x, y, z). Each element denotes the risk level of that specific voxel and is initiated to zero. The risk levels relate to the risk of a voxel being hypoxic. Higher risk levels mean higher risk of a voxel being hypoxic, while lower risk levels mean lower risk of a voxel being hypoxic. Each MS induces a risk level based on a Gaussian distribution to each voxel up to a radius of dmax. For

each MS the model loops through all voxels inside this volume and adds a risk value to the current risk level at each voxel according to the risk function

R(r, µ, σ) = exp −(r − µ) 2 2σ2 , (1)

(10)

2 GAUSSIAN RISK MODEL 2.2 Generating the risk matrix with

r =p(xM S − xvoxel)2+ (yM S − yvoxel)2+ (zM S − zvoxel)2 (2)

the Euclidean distance between the MS and the voxel. Equation 1 is a Gaussian func-tion with parameters mean µ and standard deviafunc-tion σ, and with a maximum value max(R(r, µ, σ)) = 1.

Figure 3 shows a MS that induces a risk at a voxel at distance r from the MS. The risk level at the voxel is calculated by evaluating Equation 1.

Figure 3: A MS causes a risk level at each voxel as a function of the distance r from the MS. MS: microsphere.

The two parameters µ and σ for Equation 1 are varied between each MS. The reason for this is that the brain has a complicated anatomical structure (e.g. the micro vasculature). This structure can influence how much of an effect a single lodged MS can have on the perfusion of the surrounding tissue. To represent this variability in the influence of each MS, a stochastic component is added to Equation 1. µ and σ are drawn from a Normal distribution for each MS separately, so that

µ = |N (µµ, µσ)| ·M S (3)

σ = |N (σµ, σσ)| ·M S, (4)

with M S the diameter of the MS, and N a Normal distribution

N (µµ, µσ) = 1 √ 2πµσ exp −(r − µµ) 2 2µ2 σ (5) N (σµ, σσ) = 1 √ 2πσσ exp −(r − σµ) 2 2σ2 σ , (6)

with µµand σµthe mean of the Normal distribution, and µσ and σσ the standard deviation

of the Normal distribution.

As mentioned earlier, there is a maximum distance dmax at which each voxel induces

risk. The main reason for this limitation is a limit on the computational power. Calculating and setting the risk level at each voxel is a time intensive procedure, so reducing the number of times the risk value of each voxel is updated is beneficial for the computational time.

Limiting the maximum distance has an important consequence. The risk level induced by a MS at each voxel is set according to Equation 1. This function has non-zero values below the limit r → ∞ (in the limit of r → ∞ the function goes to zero). If the limit of r = dmax, all voxels further away than r are ignored, and thus all values of

(11)

2 GAUSSIAN RISK MODEL 2.4 Statistics of output large r (since it is a Gaussian distribution), setting dmax too low can cause significant risk

values to be omitted. Setting

dmax = µµ+ 2µσ + σµ+ 2σσ (7)

ensures that at least 95 percentile of the normal distribution is taken into account. This is because for µ, 95 percentile of Equation 3 falls inside two standard deviations of the mean, and thus inside µµ+ 2µσ. The same applies to σ.

2.3 Risk matrix segmentation

Mrisk is segmented to a binary image Irisk using a threshold T , such that for every voxel i

Ii_risk= ( 0, if Mi risk< T 1, if Mi risk≥ T . (8) In Irisk the voxels with value 1 are hypoxic, and the voxels with value 0 are not hypoxic.

The aim is to find T for which the resulting HYP volume from the simulation (Vsim) equals

the HYP volume from the data (Vexp). T is calculated based on the minimisation of the

heuristic

h(T ) = (Vsim − Vexp)2. (9)

This function has exactly one minimum on the whole domain, because Vsimalways decreases

with increasing T ; when T = 0, Vsim = Vbrain, while for T > max(Mrisk), Vsim = 0. The

optimal value of T is then obtained by minimising Equation 9.

2.4 Statistics of output

Gand K functions are employed to quantify the spatial relationship between MS and HYP. Both functions are used to determine the clustering and spread of -and between- events in a spatial context (O’Sullivan and Unwin, 2010). In this study the spatial relationship between MS and HYP are examined. Therefore, two sets of events are defined: the set P containing the coordinates of all MS, and a set P containing the centroid coordinates of the HYP. By this definition, P * S.

The G function is a cumulative distribution function of the nearest-neighbour distances. The nearest neighbour distance (NND) for a point pi to the nearest event si in S is defined

as

NND(pi, S) =min(pi, si), ∀si ∈ S. (10)

Here, pi can be either part of S or a part of P (pi ∈ S or pi ∈ P).

Following from the definition of the NND in Equation 10, there are two type of G functions: homogeneous and heterogeneous. The homogeneous G function uses the NND where pi ∈ S. The heterogeneous G function, or cross-G function denoted as Gx, uses the

NND where pi ∈ P, with P * S. Since the spatial relationship between MS and HYP is

at interest in this study, and thus the relationship between P and S, only the Gx function will be treated here.

The Gx function can then be constructed summing the count of NNDs under an increasing distance d

Gx(d) = #[NND(pi, S) < d] nP

(12)

2 GAUSSIAN RISK MODEL 2.4 Statistics of output where nS and nP denote the number of events in S or P respectively.

The Gx function is used to examine the clustering of events. For example, compare two cases: events that are almost equally spaced across a domain, and events that are clustered (i.e. events are in groups, in which they are in close distance of each other). In the first case, the Gx function will have a low increase for low d. Then, halfway, it will have a rapid increase because most events have the same NND due to being equally spaced. For higher d, the function will only have a slow increase until it hits the maximum value of 1. In the second case, the Gx function will rise more rapidly for low d, because the NND is low in the event clusters. Then the Gx function will flatten for some distance, until d comes into the range of the distance between different clusters. There are many NNDs in this range for events from cluster to cluster, which means a rapid increase in the Gx function.

The K function is a cumulative distribution function of all the distances between the events, in contrast to the G function which only takes the NND in account. The function cumulatively counts the number of other events that are inside a distance d of each event. This is the same as counting all the events in a circle C(pi, d)with radius d centred around

a single event, and do this for every event, resulting in Kx(d) = PnP i=1#[S ∈ C(pi, d)] nSλ = α nPnS n X i=1 #[S ∈ C(pi, d)], ∀pi ∈ P, (12)

with α/(nPnS) a normalising factor, such that max(Kx(d)) = 1.

For the Gx and Kx functions, a squared error (SE) can be defined to quantify the difference between two functions of the same type. The SE is the squared error between two functions over their whole domain divided by the number of elements in their domain, such that SE(f(X), k(X)) = 1 n n X i=1 (f (xi) − k(xi))2, ∀xi ∈ X, (13)

with f(x) and k(x) the Gx or Kx functions of interest. X is a set of n points at which f (x) and k(x) are evaluated. The squared error (SE) for the Gx and Kx functions of the model output are then defined as

SE(Gx) ≡ SE(Gxsim, Gxexp) (14)

SE(Kx) ≡ SE(Kxsim, Kxexp), (15)

with Gxsim the Gx function of the model output and Gxexp the Gx function of the

experimental data obtained by Georgakopoulou et al. (2021), and the Kxsim the Kx

function of the model output and Kxexp the Kx function of the experimental data

obtained by Georgakopoulou et al. (2021).

The mean squared error (MSE) is defined as the average of the SE over multiple repetitions of the same simulation:

MSE(Gx) ≡ 1 k k X i=1 SE(Gx) (16) MSE(Kx) ≡ 1 k k X i=1 SE(Kx), (17)

(13)

3 METHODS & EXPERIMENTS 3.1 Implementation details

3 Methods & Experiments

3.1 Implementation details

The GRM is implemented in Python 3 (Van Rossum and Drake, 2009).2 _{For a detailed}

elaboration of the implementation details, see the README of the code included in Appendix A.

A just in time (JIT) compiler for Python 3 provided through the Numba library is used to speed up computationally intensive functions (Lam, Pitrou, and Seibert, 2015). The JIT compiler is applied as a decorator to these functions. During run-time these functions are compiled to C code by the Numba library. This can yield a significant speed-up as the Python C interpreter is bypassed by this process. This speed-up is mainly significant when used in conjunction with the NumPy library for scientific computing (Harris et al., 2020).

The brain slice for a given animal is discretized by creating a 3D risk matrix Mrisk,

whereby each element denotes a voxel in a rectangular cuboid space with coordinates (x, y, z). The units of these coordinates are voxel units. The sizes of the rectangular cuboid are set by the parameters xmax, ymax and zmax in Table 2. These parameters are given in

µm and are converted to voxel units using a scaling factor αvoxel.

αvoxel is applied to all spatial model input in order to discretize the data. All input

data has units µm and αvoxel converts these values to integer voxel units. It is important

to set αvoxelcarefully. If αvoxel is chosen to large, the model will overflow the heap memory

of Python 3 for a 8GB RAM computer. This is mainly caused by the Mrisk and Irisk,

since these two matrices contain all the voxels in the system. Since these matrices are 3D, increasing αvoxel by a factor 2 will increase the total elements in Mrisk and Irisk with

a factor 23 _{= 8}. Thus, setting α

voxel too high can quickly overflow the total available

memory.

In contrast, if αvoxel is chosen too small, there will be a significant loss of accuracy.

For example, take two coordinates x1 and x2 as input, with x1 = 100µm and x2 = 105µm.

For αvoxel = 1/6µm−1 this will result in x1 = 17voxel units and x2 = 18voxel units. For

αvoxel = 1/8µm−1 this will result in x1 = 13voxel units and x2 = 13voxel units. In the

first case there is still a difference in position for the points, but in the second case this distinction is lost and the points overlap each other.

Table 2: The constant parameters of the GRM. xmax, ymax and zmax denote the length of

the domain for the respective coordinate. αvoxel is a scaling factor applied to discretize the

domain from µm to voxel units. GRM: Gaussian risk model. Parameter Value (µm) xmax 9000

ymax 7000

zmax 500

αvoxel 1₆

The MS positions and sizes are stored in a matrix of size (#MS, 4), with the first column denoting the diameters of the MS and the other three columns the positions in (x, y, z) coordinates, all in voxel units.

(14)

3 METHODS & EXPERIMENTS 3.2 Experiments The HYP centroids are stored in a matrix of size (#HYP, 3), with the columns denoting the positions in (x, y, z) coordinates in voxel units.

As explained in subsection 2.3, Equation 9 is minimized to find the optimal value of T to segment the Mrisk to the Irisk. This minimisation problem is solved using

SciPy’s optimize.minimize_scalar function (Virtanen et al., 2020). This function finds the local minimum of a scalar function. Since Equation 9 has exactly one minimum, optimize.minimize_scalar is guaranteed to find the optimal value of T .

The Irisk contains the hypoxic regions calculated by the model. The properties of these

regions are calculated using scikit-image’s measure.label and measure.regionprops functions (van der Walt et al., 2014). The last function returns a props objects, which contains information about the number of HYP, the total HYP volume, and the coordinates of the centroids of the HYP.

3.2 Experiments

The GRM has four variable parameters: µµ, µσ, σµ and σσ. The aim is to find the

parameter set for which the spatial model output fits the HYP data of each animal the best. Since the input per animal differs significantly (see Table 1), the parameters are optimised for each animal individually. The statistic to judge the model outcome to the data is the Gx function. Therefore, the objective is to find the set of parameters for which MSE(Gx) is minimal, as defined in Equation 16.

Since the system is complex and may have many local minima for MSE(Gx), a stochastic minimisation algorithm is used: SciPy’s optimize.dual_annealing with no_-local_search = True. This is a simulated annealing algorithm that aims to find the global minimum of a given function, in this case the MSE(Gx). The algorithm iterates for a given number of steps. In each step of optimize.dual_annealing the simulation is performed k = 10 times to calculate the MSE(G). The maximum number of iterations is set to maxiter = 100. The algorithm returns the optimised parameters for which the MSE(G) is minimal. Since simulated annealing is a stochastic global minimisation algorithm, there is no guarantee that the global minimum will be found. That is why it is important to set maxiter high enough to increase the chance of finding the global minimum. However, even if the global minimum is not found, after a limited amount of iterations, the solution can be assumed to be relatively close to the global minimum.

This minimisation procedure results in four unique parameter values for each animal. To validate the universality of the optimised parameter sets (i.e. how well the parameter set of one animal works to simulate the HYP for the different animals), the optimised parameters are cross-applied to each animal. This means that the optimised parameters for animal 1 are applied to animals 1 to 6, the optimised parameters for animal 2 are applied to animals 1 to 6, et cetera. The output is judged on the spatial metrics MSE(Gx), σ(MSE(Gx)), MSE(Kx), and σ(MSE(Kx)). σ(MSE(Gx)) and σ(MSE(Kx)) are the standard deviations over k simulations of MSE(Gx) and MSE(Kx) respectively.

A one-at-a-time (OAT) sensitivity analysis is performed on the output parameters to test the sensitivity of each parameter. This can give insight on the influence of single parameters. During OAT sensitivity analysis a single parameter is varied, while keeping the other parameters at a (constant) baseline. The output of this analysis is again judged on the MSE(Gx), σ(MSE(Gx)), MSE(Kx), and σ(MSE(Kx)).

The optimal threshold for each animal for their optimised parameters averaged over 100 repetitions is examined to decide if the segmentation based on fitting the threshold to

(15)

3 METHODS & EXPERIMENTS 3.2 Experiments the total volume of HYP is necessary. If the threshold varies significantly between animals, it means that the segmentation method was necessary to obtain valid results for the other metrics. If the threshold does not vary significantly, it means that one universal threshold can be introduced in the model, removing the necessity of knowing the total HYP volume beforehand.

(16)

4 RESULTS 4.1 Parameter optimisation

4 Results

Figure 4: A MIP of two MS (white dots), who both causes a risk area according to a Gaussian distribution. The risk levels are coloured in red, while black denotes areas free of risk. MIP: maximum intensity picture, MS: microsphere.

Figure 4 shows an example of two MS from Figure 2 causing a risk area around them according to the risk function in Equation 1. Both MS have a different value for µ and for σ, drawn from Equation 3 and Equation 4 respectively. This can be observed by the fact that the upper MS has a more shallow "ring" of high risk levels than the lower MS. The thickness of this "ring" is governed by σ, the width of the Gaussian distribution.

The MS have an overlap of their risk levels in the region directly between the MS. These risk levels are added together, forming a total risk level at each voxel. From Figure 4 it can be observed that the highest risk levels are found in this region. As mentioned earlier, a higher risk level in an area means more chances of having a hypoxic region in that areas. This is an example of how MS can work together to form HYP.

4.1 Parameter optimisation

The optimised sets of parameters obtained through minimising the MSE(Gx) using a simulated annealing algorithm for all animals are presented in Figure 5. µσ and σσ both

have low values for most of the animals, which means that there is a low stochastic component of the model, since these two parameters govern the stochastic component of the model. µµ is higher than σµ, except for animal 6. The optimised parameter set for

animal 6 also has a relatively high value for σσ. There is a large variation of the parameter

values between animals.

The average number of HYP generated by the model for the optimised parameter set of each animal is presented in Figure 6. Strikingly, the number of HYP for animal 1 has very high value compared to the number of HYP from the experimental data. The GRM predicts the number of HYP better for the five other animals, although animal 6 also deviates significantly from the data. In general, there is significant variation in the

(17)

4 RESULTS 4.1 Parameter optimisation

Figure 5: The optimised parameters of the GRM for each animal. There is a significant difference in the optimised parameters for each animal. µµ is higher than σµ except for

animal 6. Both standard deviations µσ and σσ have relatively low values compared to the

means µµand σµ. GRM: Gaussian risk model.

Figure 6: The average number of HYP using the optimised parameters for each animal. The error bars show the standard deviation over 100 repetitions. HYP: hypoxia.

(18)

4 RESULTS 4.1 Parameter optimisation number of HYP generated by the model, except for animals 4 and 5. This indicates that the number of HYP is sensitive under the stochastic component of the model: the µσ and

σσ parameters.

Figure 7: The optimised parameter sets of each animal (on the vertical axes) are used to perform the simulation for all animals (on the horizontal axes). Moving clockwise and starting from the left top corner, the MSE(Gx), σ(MSE(Gx)), MSE(Kx), and σ(MSE(Kx)) are plotted. Each simulation is repeated 100 times. MSE: mean squared error, σ: standard deviation.

Figure 7 contains the results of the simulation of the data of each animal, each performed for all six optimised parameter sets. The four resulting matrices show how well the spatial functions of the simulation fit the experimental data. Out of all the parameter sets, the optimised parameters of animal 1 fits the MSE(Gx), σ(MSE(Gx)), MSE(Kx), and σ(MSE(Kx)) best of all animals. At the same time, individual animals can have a better fitted parameter set (e.g. animal 5 shows the best fit for its own optimised parameter set).

Moreover, Figure 7 shows that each animal does not always have the best fit for MSE(Gx), σ(MSE(Gx)), MSE(Kx), or σ(MSE(Kx)) with its own optimised parameters. For animal 3, 4 and 6 the MSE(Gx) is even the highest for their own optimised parameters than any of the other parameter sets. This is explained by the simulated annealing algorithm used to find these parameters. As mentioned earlier, the algorithm is not guaranteed to find the global minimum in a limited number of iterations. Since the global minima for at least 3 of the animals were not found, not enough iterations of the simulated annealing algorithm were used.

Finally, the optimised parameter set of animal 3 has poor performance for all four metrics MSE(Gx), σ(MSE(Gx)), MSE(Kx), and σ(MSE(Kx)) compared to the parameter sets of the other animals. This might be explained by the fact that animal 3 has the lowest amount of HYP and the lowest hypoxic volume in the experimental data (Table 1). This differs significantly with the other animals, which have both higher numbers of HYP and

(19)

4 RESULTS 4.2 Spatial analysis higher HYP volume.

4.2 Spatial analysis

The Gxsim for each animal, with their respective optimised parameters, are plotted with a

standard deviation over 30 repetitions against the Gxexp functions of the data in Figure 8.

These plots show significant differences between how well the model output of each animals fits the data.

Especially the Gxsim of animals 3 and 6 show a poor fit to Gxexp as seen in Figure 8.

The Gxexp for these animals has a steeper increase of the cumulative fraction than Gxsim,

meaning that many nearest-neighbour HYP are closer to their MS than in the model output. This also indicates that in reality there is less variation in the NND between MS and HYP than in the model output, since in the former there is less spread of the Gxexp

over the NND than for the Gxsim in the latter.

The standard deviation of the Gxsim is generally low, but animals 3 and 4 have notable

variance, of which animal 4 has the largest variance. Animal 4 has the highest µσ parameter

of all animals (Figure 5). This parameter potentially has a large influence on the NND of a MS to a HYP, because it governs the variance of the distance of the maximum of the risk level caused by a MS. A MS is expected to have a significant influence on the closest HYP, because the MS induces a risk area closely around itself. A higher µσ then leads to

more randomness in the NND of MS to HYP. Thus, the Gxsim, a function of the NND,

varies more for higher µσ.

Figure 8: The Gx functions of the NND between MS and HYP for the model output and the experimental data for each of the six animals (counted from left to right, top to bottom). The standard deviation is shown for the model output over 30 repetitions. HYP: hypoxia, MS: microsphere, NND: nearest neighbour distance.

The Kxsim for each animal, with their respective optimised parameters, are plotted

with a standard deviation over 30 repetitions against the Kxexp of the experimental data in

(20)

4 RESULTS 4.3 Sensitivity analysis in general between 2000 µm and 5000 µm distance between MS and HYP. For all animals, except animal 4, the Kxsim is slightly skewed to the left in this interval compared to the

Kxexp. This means that the HYP in the model output are slightly more clustered than

the HYP in the data.

It is noteworthy to mention that animal 4 has the highest standard deviation of all the animals for Kxsim in Figure 9, as it also did for Gxsim in Figure 8. This can again

be explained by the large values for µσ and the resulting higher variance in the model

output. Interestingly, the Kxsim of animal 4 seems to have a different shape of the graph

than Kxexp. The characteristic curve of the Kxexp, where at higher distances the rate of

increase of the cumulative fraction reduces, does not seem to appear in the Kxsim. Instead,

a more straight linear curve is observed.

Figure 9: The Kx functions of the distances between MS and HYP for the model output and the data for each of the six animals (counted from left to right, top to bottom). The standard deviation is shown for the model output over 30 repetitions. HYP: hypoxia, MS: microsphere.

4.3 Sensitivity analysis

This section treats the results of the OAT sensitivity analysis for each of the four model parameters. Each time during this analysis three parameters are fixed at the optimised values obtained in Figure 5 and the remaining parameter is varied over a range of values. The output of this analysis is judged on the differences between the Gx and Kx functions of the model output and experimental data: MSE(Gx) and MSE(Kx), and their standard deviations σ(MSE(Gx)) and σ(MSE(Kx)) respectively.

First, Figure 10 shows the influence of the µµparameter on the spatial metrics. Animals

1, 2, 4 and 5 show clearly that low values for µµ leads to a worse fit to the MSE(Gx)

compared to higher values for µµ. There is, however, no clear value for which the parameter

(21)

4 RESULTS 4.3 Sensitivity analysis

Figure 10: Sensitivity of µµ for the optimised parameter set for each animal. Moving

clockwise and starting from the left top corner, the MSE(Gx), σ(MSE(Gx)), MSE(Kx), and σ(MSE(Kx)) are plotted. The standard deviations are calculated using 10 repetitions. MSE: mean squared error, σ: standard deviation.

Animal 3 shows less variation of MSE(Gx) over the whole range of µµ values than the

other animals, although µµ= 0, which is the lowest possible value, still performs worse

than any higher value of µµ.

The MSE(Gx) of animal 6 shows the opposite trend compared to the other animals, whereby lower µµ values reduce the MSE(Gx) instead of increase the MSE(Gx). This

means that the optimal µµ for animal 6 is lower than for the other animals according to

the MSE(Gx). This is in accordance with the optimised parameters of animals 6, where µµ has a very low value compared to the other animals (Figure 5).

MSE(Kx) does not show any clear relationship with µµ. The values of MSE(Kx) are

generally lower than for MSE(Gx), meaning that µµ mostly influences how well the Gx

functions of the model output and the experimental data fit, and has less influence on how well the Kx functions of the model output and the experimental data fit. This indicates that µµ has mainly a local spatial influence, since the Gx function describes the NND

between MS and HYP, opposed to the Kx describing the distances of each MS to all HYP. The Gx function is thus more influenced by local events, while Kx is less sensitive to local events.

Both σ(MSE(Gx)) and σ(MSE(Kx)) do not show any clear relationship between µµ and the respective metric. There is slightly more variation of σ(MSE(Gx)) than of

σ(MSE(Kx)), but both have low values compared to MSE(Gx). That both metrics do not have a clear relationship with µµ is expected, since varying µµ does not change the

stochasticity of the model, meaning that µµ does not influence the variance of the model

output.

Second, Figure 11 shows the influence of the µσ parameter on the spatial metrics.

MSE(Gx) shows a trend whereby it increases for increasing µσ, meaning that the µσ

(22)

4 RESULTS 4.3 Sensitivity analysis parameter is 0 µm. Since µ is drawn from Equation 3 where the standard deviation equals µσ, setting µσ = 0 means that Equation 3 changes to a deterministic formula. Thus, there

is no variance in the value of µ between different MS for µσ = 0µm.

Figure 11: Sensitivity of µσ for the optimised parameter set for each animal. Moving

The MSE(Kx) as a function of µσ shows only a small effect for animal 4, where

decreasing µσ leads to a decrease in MSE(Kx). For the other animals, there is no

significant difference caused by µσ. This is in line with the MSE(Kx) in Figure 10, where

the MSE(Kx) seemed to be mostly independent of µσ compared to the MSE(Gx). This

can again be explained by the same argument advanced earlier, that the µσ parameter has

more influence locally than over larger distances.

σ(MSE(Gx)) shows a small negative correlation with µσ. Smaller µσ seem to decrease

σ(MSE(Gx)) and thus lead to smaller variance of MSE(Gx). There is one outlier for animal 4 around µσ = 44µm.

Similar to MSE(Kx), σ(MSE(Kx)) as a function of µσ shows only a small effect for

animal 4, where decreasing µσ leads to a decrease in σ(MSE(Kx)).

Next, Figure 12 shows the influence of the σµ parameter on the spatial metrics. A

few interesting observations can be made about the relation of MSE(Gx) with σµ. The

first is that animals 1, 2, 4 and 5 have a very high MSE(Gx) for σµ, around 0 µm to

30µm. Since σµ has low values here, the width of the risk function in Equation 1 is small,

meaning that each MS only causes a very small and thin risk area. Since the total HYP volume of the model output is fitted to the experimental data, most of these areas are necessarily denoted as HYP, leading to an inaccurate Gxsim and thus a large MSE(Gx).

The second observation is that increasing σµ to high values can also lead to a higher

MSE(Gx). Making σµ very high effectively causes the risk function in Equation 1 to have

widely distributed values. Since for high values of σµ, µµ is close to or smaller than σµ, the

(23)

4 RESULTS 4.4 Threshold

Figure 12: Sensitivity of σµ for the optimised parameter set for each animal. Moving

chance of a HYP to be very close to the MS, reducing the NND between MS and HYP. Since the Gx function is a function of the NND between MS and HYP, the MSE(Gx) increases, precisely as observed in Figure 12.

Interestingly, σ(MSE(Gx)) has a few very high values around σµ = 20µm − 40 µm

compared to other values of µσ.

The MSE(Kx) and σ(MSE(Kx)) both do not show any clear relationship with σµ.

The largest variation in their values is (again) found for animal 4.

Finally, Figure 13 shows the influence of the σσ parameter on the spatial metrics. The

exact relationship of MSE(Gx) and σσ varies between the animals. However, in general

there seems to be a small trend where MSE(Gx) has lower values for smaller σσ. Only

animal 6 seems to have an opposite trend.

σ(MSE(Gx)) does not show a clear trend as a function of σσ, but most animals,

especially animal 4, have several significant values for higher values of σ(MSE(Gx)). This is expected, since σσ describes the variability of the σ parameter from Equation 4 for

each MS. Having a high value for σσ, can introduce more variance between subsequent

simulations.

MSE(Kx) and σ(MSE(Kx)) both do have low values and vary little as a function of σσ compared to MSE(Gx) for all animals, except animal 4. For this animal an increase of

σσ results in an increase of MSE(Kx) and σ(MSE(Kx)).

4.4 Threshold

The optimal threshold for segmenting the HYP from the risk matrix, found using the heuristic in Equation 9, varies significantly between animals. The highest threshold is more than twice as high as the lowest threshold in Figure 14. An interesting notion here is that

(24)

Figure 13: Sensitivity of σσ for the optimised parameter set for each animal. Moving

all the thresholds are higher than 1. Recall, following from the definition of Equation 1, that max(R(r, µ, σ)) = 1. This means that a single MS can only induce a maximum risk level of 1 to any voxel. Since the threshold is defined as the minimum risk level that is required for a voxel to be hypoxic, a single MS cannot cause HYP following from this definition. For animal 1, 2, 3 and 6 at least 2 MS are needed to create a HYP, and for animal 4 and 5 at least 3 MS are needed to create a HYP.

The total HYP volume in the model is a continuous function of the threshold and decreases over the whole domain. This decrease is, however, not linear. In Figure 15 there are multiple tipping points visible for each animal. These tipping points generally seem to occur around multiples of one. This can, again, be explained by the maximum risk value induced by a single MS: max(R(r, µ, σ)) = 1. If the threshold passes the value of 1, suddenly a single MS cannot increase a HYP anymore, which leads to a sharp decrease in the potential total HYP volume. This transition again occurs when the threshold passes the value of 2, when two MS cannot cause a HYP anymore without help of more MS. This trend is, to a lesser extend, also observed around threshold 3, 4 and 5.

Another important feature from Figure 15 is that the maximum risk value differs greatly between the animals. This happens because the maximum threshold value is equal to the highest risk level in Mrisk. For example, animal 5 has almost 2.5 times higher

maximum risk level than animal 4. This shows that there are large differences between the model output for the different animals.

(25)

Figure 14: The average threshold for the optimised parameters for each animal. Observe that all the thresholds have value higher than 1, meaning that in order to create a HYP, at least two MS need to cooperate. The error bars show the standard deviation over 100 repetitions. MS: microsphere.

Figure 15: The threshold as a function of the total HYP volume for each animal using their optimised parameters. Around the integer values of the thresholds, transitions of the graphs are visible. These values correspond with the maximum risk level of 1 caused by a single MS. MS: microsphere, HYP: hypoxia.

(26)

5 DISCUSSION

5 Discussion

In this study, a new stochastic model was presented to simulate HYP caused by MS in a rat model of the brain. The model was applied to the data of MS and HYP of six rats from the experiments of Georgakopoulou et al. (2021).

The model output showed significant variability for all metrics between the different animals. This reflects the differences in the input data of each animal. There is more than a four-fold difference between the lowest and highest number of MS across the animals. This can significantly alter the MS density (Georgakopoulou et al., under review). Because the GRM assumes an additive linear relationship between MS, decreasing the MS density should lead to less cooperation between MS since the MS are more spread out compared to higher densities. In the same fashion, increasing the MS density should lead to more cooperation between MS.

This hypothesised positive relationship between the number of MS and cooperation between MS is, however, not completely supported by the GRM. The total number of HYP generated by the model for each animal does imply that this relation exists. The MS of animal 3, with the lowest MS count, generated the least number of HYP. The other animals, with higher MS count generated more HYP. However, the optimal threshold for each animal to segment the Mrisk to the Irisk, showed that animal 3 had a similar

threshold as animals 1 and 6 and even higher than animal 2. These three animals all have significantly more MS than animal 3. As mentioned earlier, the threshold gives a lower limit on the number of MS that are needed to create a HYP. For all of the animals 1, 2, 3 and 6 at least two interacting MS are needed to create a HYP. This implies that the cooperation between MS does not change with the total number of MS.

The optimised parameters that were obtained for each animal using simulated annealing do not reflect the optimised parameter set for each animal. However, some of the optimised parameter sets translate reasonably well between animals. Especially the optimised parameter set of animal 1 gives a low error in the spatial functions when applied to the simulation of the other animals. These results confirmed that it is not guaranteed that the simulated annealing algorithm finds the optimal parameters in the limited number of iterations performed.

Interestingly, there does not seem to be a direct connection between the standard deviation of the number of HYP and the standard deviation of the Gx function. Recall that animal 1 has significant variability in the number of HYP generated over multiple simulation runs. Take the case of a simulation that has a low number of HYP for this animal. Then, if the Gx function between MS and HYP is constructed, there are a (relatively) low number of HYP available as nearest-neighbour to each MS. Take now the case that a simulation results in a high number of HYP. Suddenly, there are more HYP candidates as nearest-neighbour to each MS. Since the resulting Gx functions do not vary much between these two cases, it means that the NNDs between MS and HYP do not differ much either. This implies that the HYP in the first case gets broken up in smaller HYP in the second case (the total HYP volume stays equal, thus more HYP means smaller HYP).

The fact that the HYP appear to get broken up in smaller pieces without affecting the Gx function, also indicates that the number of HYP is not a very robust metric. For example, imagine two small HYP with each the size of a few voxels, separated from the other by just one voxel. Is it justifiable to treat them as two different HYP, or is this separation just caused by numerical errors and mathematical approximations? It might

(27)

5 DISCUSSION 5.1 Limitations & Future recommendations even be, that if the threshold would have been just slightly higher, suddenly the former two HYP would appear to be one. This issue also comes into play when treating the segmentation of HYP from brain images. It is not always clear if, and where, there is a division between different closely situated HYP.

The sensitivity analysis shows that for a majority of parameter settings MSE(Gx) > MSE(Kx). Also, MSE(Kx) has less variability between different parameter settings than MSE(Gx). Since the Gx function looks at the NND between MS and HYP, changing the parameters directly impact the Gx function. This is because a MS creates a risk areas locally, which means that it has a large influence on where the closest HYP is created. Over larger distances however, the influence on the risk level of that area caused by this MS reduces. Other MS might be located in these areas which exert a larger influence on the HYP created there. Changes in the simulation parameters will in this case not significantly alter the distance from the MS to further HYP. The Kx function captures this latter part; the distance from the MS to all HYP, of which the majority are relatively far away from the MS. Thus the Kx function is expected to vary less than the Gx function, which is indeed observed.

The sensitivity analysis also indicates that the model is more sensitive to the µµ and

σµ parameters than to the µσ and σσ parameters. The last two parameters seem to have a

small negative correlation with the spatial metrics, meaning that setting them to non-zero values generally leads to worse model performance. Therefore, it is possible to remove them from the model by setting their values to zero. This will mean that Equation 3 and Equation 4 do not depend on a normal distribution anymore, but instead will be deterministic formulas proportionate to µµ and σµ respectively. This has an important

consequence, because µσ and σσ represent the stochastic component of the model. By

removing these parameters from the model, the model becomes deterministic.

Removing µσ and σσ has as consequence that the variability between different MS is

taken out of the model. The assumption was made in subsection 2.2 that each MS induces a slightly different risk area, since the underlying anatomical structure varies between the positions where the MS are lodged. If the formulas of µ and σ are changed from stochastic to deterministic ones, every MS will induce similar risk areas.

The average threshold differs significantly between each animal, as does the maximum risk level in the Mrisk. This indicates that there are significant differences between

simulations of the HYP of the different animals. The threshold as a function of the total HYP volume showed that there is a minimum number of two or three MS needed to induce HYP in the model. The question is whether this conclusion can be translated to the in-vivo rat data of Georgakopoulou et al. (2021). It can be argued that, since the spatial relationships between MS and HYP are similar for both, that if the MS work together to create HYP in the model, this implies that the same holds true for the in-vivo system.

5.1 Limitations & Future recommendations

A few limitations of the model need to be mentioned. First of all, the model uses a discretized matrix to represent the brain slices. This causes a loss in accuracy compared to a continuous model. During discretization the brain slice is divided into voxels, with each voxel a length of αvoxel times the length of the corresponding brain volume in µm.

The data is supplied with an accuracy smaller than µm, thus there is an error introduced by the discretization.

(28)

5 DISCUSSION 5.1 Limitations & Future recommendations cuboid, opposed to the real shape of a brain. Since HYP can be created in the whole matrix in this model, HYP can be (partially) situated outside the real brain boundary. To compensate for this, a mask can be applied that removes these HYP which are situated outside of the brain. However, data on the exact brain boundaries are needed to implement this.

Second, the data describes the MS and HYP in a 500 µm slice of the rat brains. However, there also exist MS and HYP outside of this slice, in other parts of the brain. The fact that there are MS outside this brain slice, which are not taken into account in this model, is important. These MS can still cause HYP that lay inside the brain slice. By not taking this into account, a higher emphasis is placed on the role of the MS located in the brain slice. To counter this, a balancing mechanism should be introduced. For example, periodic boundaries at the cut sides of the brain slice may create a more realistic contribution of MS to the system. An alternative is to create statistically accurate MS in-silico over a full brain volume of the mouse by extrapolating the data.

Furthermore, the total HYP volume is supplied to the model in the data input. This information is needed to segment the Mrisk to the Irisk. However, ideally any information

about the HYP would be unnecessary before running the model. It is preferred to perform the simulation based on just the data of the MS alone. This projects the physical systems better, since in reality the MS cause the HYP, leading to a total HYP volume. It does not happen the other way around i.e. the MS and total HYP volume cause HYP.

Doing the segmentation without any information on the total HYP volume is quite complicated though, as demonstrated in subsection 4.4. The optimal threshold varies significantly between animals, as does the maximum risk level of the Mrisk, as demonstrated

in Figure 15.

Finally, the model assumes that there is an linear additive relationship between the risk levels of the MS. The risk levels of different MS at a voxel are added together linearly, but this is not necessarily a true representation of reality. The micro vasculature, where the MS are lodged, has a complex layout, where the vessels can intersect other vessels and form loops. This means that having more MS lodged in the vessel bed, does not per definition increases the risk linearly. It could for example also increase exponentially, or logarithmically.

The micro vasculature can be taken into account in various ways. One way is to model the micro vasculature explicitly i.e. the vessels themselves, their connectivity and blood flow. This is possible to do, and actually applied to micro embolisms by Xue et al. (2021). Another possibility is to look at the micro vasculature density and its spatial relationship with MS and HYP.

(29)

6 CONCLUSION

6 Conclusion

The spatial relationship between MS and HYP generated by the model has similar characteristics as the experimental data, but there is no complete overlap between their spatial functions. The maximum risk levels induced by the MS differs significantly between the animals. It is concluded, from the threshold which determines which brain area is hypoxic or not, that multiple MS need to induce a high enough risk level to the same area, in order for that area to become hypoxic. This poses the question whether this is also the case for the causation of HYP by micro-emboli in human patients.

Although caution should be taken in translating the results from an in-silico model of rodents to humans (because of obvious anatomical differences), the results of this study suggests that it is necessary to design a stent retriever for the EVT procedure that reduces the fragmentation of the emboli into smaller pieces as much as possible. Having a lower amount of micro-emboli lodged in the brain of patients can have a positive impact on their clinical outcome, since multiple MS are needed to create HYP.

(30)

6 REFERENCES

References

Bhaskar, S., Stanwell, P., Cordato, D., Attia, J., & Levi, C. (2018). Reperfusion therapy in acute ischemic stroke: dawn of a new era? BMC neurology, 18 (1), 1–26. https: //doi.org/10.1186/s12883-017-1007-y

Chueh, J.-Y., Puri, A. S., Wakhloo, A. K., & Gounis, M. J. (2016). Risk of distal emboliza-tion with stent retriever thrombectomy and ADAPT. Journal of neurointervenemboliza-tional surgery, 8 (2), 197–202. https://doi.org/10.1136/neurintsurg-2014-011491

Donkor, E. S. (2018). Stroke in the century: a snapshot of the burden, epidemiology, and quality of life. Stroke research and treatment, 2018. https://doi.org/10.1155/2018/ 3238165

Georgakopoulou, T., van der Wijk, A.-E., Bakker, E. N., & van Bavel, E. (2021). Recovery of Hypoxic Regions in a Rat Model of Microembolism. Journal of Stroke and Cere-brovascular Diseases, 30 (6), 105739. https://doi.org/10.1016/j.jstrokecerebrovasdis. 2021.105739

Georgakopoulou, T., van der Wijk, A.-E., Bakker, E. N., & van Bavel, E. (under review). Quantitative 3D analysis of tissue damage in a rat model of microembolization. Goyal, M., Demchuk, A. M., Menon, B. K., Eesa, M., Rempel, J. L., Thornton, J., Roy, D.,

Jovin, T. G., Willinsky, R. A., Sapkota, B. L., Dowlatshahi, D., Frei, D. F., Kamal, N. R., Montanera, W. J., Poppe, A. Y., Ryckborst, K. J., Silver, F. L., Shuaib, A., Tampieri, D., . . . Hill, M. D. (2015). Randomized Assessment of Rapid Endovascular Treatment of Ischemic Stroke [PMID: 25671798]. New England Journal of Medicine, 372(11), 1019–1030. https://doi.org/10.1056/NEJMoa1414905

Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del R’ıo, J. F., Wiebe, M., Peterson, P., . . . Oliphant, T. E. (2020). Array programming with NumPy. Nature, 585 (7825), 357–362. https://doi.org/10.1038/s41586-020-2649-2

INSIST. (2021). About us. Retrieved May 16, 2021, from https://www.insist-h2020.eu/ images/sliderimages/plaatje-3.jpg

Konduri, P. R., Marquering, H. A., van Bavel, E. E., Hoekstra, A., Majoie, C. B. L. M., &. (2020). In-Silico Trials for Treatment of Acute Ischemic Stroke. Frontiers in Neurology, 11, 1062. https://doi.org/10.3389/fneur.2020.558125

Lam, S. K., Pitrou, A., & Seibert, S. (2015). Numba: A llvm-based python jit compiler. Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, 1–6. https://doi.org/10.1145/2833157.2833162

O’Sullivan, D., & Unwin, D. J. (2010). Geographic information analysis and spatial data. Geographic information analysis, 130–137. https://doi.org/https://doi.org/10.1002/

9780470549094

Phipps, M. S., & Cronin, C. A. (2020). Management of acute ischemic stroke. bmj, 368. https://doi.org/10.1136/bmj.l6983

Samuel, J., & Franklin, C. (2008). Hypoxemia and hypoxia. Common surgical diseases (Pages 391–394). Springer. https://doi.org/10.1007/978-0-387-75246-4_97

van der Walt, S., Schönberger, J. L., Nunez-Iglesias, J., Boulogne, F., Warner, J. D., Yager, N., Gouillart, E., Yu, T., & the scikit-image contributors. (2014). scikit-image: image processing in Python. PeerJ, 2, e453. https://doi.org/10.7717/peerj.453 Van Rossum, G., & Drake, F. L. (2009). Python 3 Reference Manual. CreateSpace.

(31)

6 REFERENCES

van Veluw, S. J., Shih, A. Y., Smith, E. E., Chen, C., Schneider, J. A., Wardlaw, J. M., Greenberg, S. M., & Biessels, G. J. (2017). Detection, risk factors, and functional consequences of cerebral microinfarcts. The Lancet Neurology, 16 (9), 730–740. https://doi.org/10.1016/S1474-4422(17)30196-5

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., . . . SciPy 1.0 Contributors. (2020). SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17, 261–272. https://doi.org/10.1038/s41592-019-0686-2

Wardlaw, J. M., Murray, V., Berge, E., & Del Zoppo, G. J. (2014). Thrombolysis for acute ischaemic stroke. Cochrane database of systematic reviews, (7). https://doi.org/10. 1002/14651858.CD000213.pub3

World Health Organization. (2020, December 9). The top 10 causes of death. Retrieved May 10, 2021, from https://www.who.int/news-room/fact-sheets/detail/the-top-10-causes-of-death

Xue, Y., El-Bouri, W. K., Józsa, T. I., & Payne, S. J. (2021). Modelling the effects of cerebral microthrombi on tissue oxygenation and cell death. bioRxiv. https: //doi.org/10.1101/2021.01.16.426717

(32)

A README A.2 Usage

A

README

This is a description of the code for the Gaussian risk model. This model creates hypoxia (HYP) based on data of microspheres (MS, or ms in the code) of rats.

In total the data of 6 rats are used with the identifiers: 1828, 1830, 1849, 1858, 1869 and 1872.

All input data has micrometer units. Units are converted to voxel units by dividing them by a scaling factor voxel_size and rounding them to the nearest integer.

COMMON ABREVIATIONS & NOMENCLATURE: • GRM: Gaussian risk model

• HYP: hypoxia

• MS / ms: microsphere(s)

• NND: nearest-neighbour distance • np: NumPy

• risk: commonly used in the code as a substitute for HYP

NOTE: This README gives a general introduction to each source file, and summarises the important functions. For explicit coding details and elaborate function descriptions view the comments and doc-strings of the functions in the source files.

A.1 Author

Sjoerd Terpstra

sjoerd.terpstra@student.uva.nl

A.2 Usage

Python 3 is required to run the model. Python version >=3.8 is recommended, older versions might work, but that is not guaranteed.

A.2.1 Program setup

These are the setup instructions for Linux and MacOS: • Install pip3: sudo apt-get install python3-pip

• Install virtualenv using pip3: sudo pip3 install virtualenv • Create virtualenv: python3 -m venv myenv

• Active the virtual environment: source myenv/bin/activate • Install required packages: pip3 install -r requirements.txt • Run model: python3 main.py

• Deactivate the virtual environment after usage: deactivate

For Windows make sure to install all the required packages in requirements.txt. It is recommended to use anaconda for this, especially for correct installation of the Numba library (although pip3 should also work on Windows).

(33)

A README A.3 Modules

A.3 Modules

The model is run from main.py, which is the most important file. In this file all the simulations and model runs are executed. What follows is a short description per module, and their most important functions/functionalities.

A.3.1 data_merge_script.py

Used as a stand-alone script to merge the data sets of the MS coordinates (data/Coordinates_-microspheres/) and the corresponding diameters of each MS (data/Diameters_micro-spheres/) to new data sets containing the information about the MS coordinates and diameters in data/Data_microspheres/. These last data sets are used to load any data of the MS in the simulation.

A.3.2 hypo.py

This file contains the function experimental_hypo_centroids() to load the centroids of the experimentally obtained HYP. It reads the files located in data/Coordinates_-hypoxic_centroids/. Each file contains the centroid coordinates in micrometers of each animal in .csv files. The centroids are stored in a np array of size (#HYP, 3), with the 3 columns denoting the x, y and z coordinate of the centroid. These coordinates are converted to voxel units (see parameters.py).

A.3.3 main.py

This is the main file, from which the program is run. All other files, except data_merge_-script.py, are called from here. This file contains all the different simulations (e.g. parameter optimisation, generating the threshold curve et cetera) in separate function. Each of these functions is documented in the code. Only a few functions are highlighted here.

• main(): this is the entry point of the program, and the function that is execute if the file is executed through python3 main.py.

• run_gaussian_risk_model(): this function executes the GRM itself. It takes two sets of parameters:

– pars_model: this class object contains all the contant model parameters, as defined in parameters.py.

– pars_simul: this is a dictionary containing the following parameters:

∗ animal: index of animal, corresponding to the entry in pars_model.animals, for which the simulation is executed.

∗ mu_mean (see risk.py). ∗ mu_stddev (see risk.py). ∗ sigma_mean (see risk.py). ∗ sigma_stddev (see risk.py).

∗ threshold_hardcoded (see segment.py). ∗ threshold_multiplier (see segment.py).

(34)

A README A.3 Modules It consists of three basic steps:

1. Load the ms data using ms.load_ms().

2. Generate the risk matrix by calling risk.risk_gaussian().

3. Segment the risk matrix to a risk image containing the hypoxic regions by calling segment.segment_risk_matrix().

The function returns the loaded ms, the risk matrix, and the risk image. Most simulations call this function, since this forms the core of the model.

• load_pars_simul(): this function loads the pars_simul saved by the optimisation algorithm in the results/ directory. Each animal has its own optimised parameter set.

A.3.4 ms.py

This file contains the function load_ms() to load the diameters and coordinates of all the MS for a single animal. This function has an argument ms_type, with which the type of MS that is loaded can be specified. These two types of MS can be loaded:

• Random MS: these are uniformly random generated MS with a specified diameter in parameters.py. This function does not take into account the real brain boundaries, only the rectangular cuboid with sizes specified in parameters.py. It is recommended to not use this option in the current form.

• Data MS: these are the MS from the 6 rats. The diameters and coordinates are stored in .csv files located in data/Data_microspheres/. These values are converted to voxel units (see parameters.py).

The MS are stored in a np array of size (#MS, 4), with the first column denoting the diameters and the other 3 columns the x, y and z coordinates. Only the MS from a single specified animal are returned by load_ms().

A.3.5 parameters.py

Contains all the constant parameters for the whole model in a single class Parameters. This is invoked only once in the top of the main() function in main.py and saved in the object pars_model. The most important parameters used are listed here:

• animals: list of the identifiers for the 6 rats.

• hypo_volume_tot: list of the total HYP volume for the 6 rats. The indices correspond with the associated animal in animals.

• n_hypo_tot: list of the total number of HYP for the 6 rats. The indices correspond with the associated animal in animals.

• path_data_ms: path to files where the MS data is stored. • path_data_hypo: path to files where the HYP data is stored.

• voxel_size: scaling factor to convert micrometers to voxel units. This is the inverse of the parameter αvoxel mentioned in the thesis.

• brain_size_micron: the length of the x, y, z axis of the rectangular cuboid around the slice of the rat brain that is simulated. Units are in micrometers.

• brain_size_voxel: the length in the x, y, z axis of the rectangular cuboid around the slice of the rat brain that is simulated. Units are in voxel units.

• ms_size_micron: the three different possible MS diameters in micrometers. • ms_size_voxel: the three different possible MS diameters in voxel units.

In-silico model of hypoxic areas caused by micro-embolisms in the rat brain