Adapting submodules of the in-stent restenosis model to enable uncertainty quantification using MUSCLE 3

(1)

Bachelor Informatica

Adapting submodules of the

in-stent restenosis model to

en-able uncertainty quantification

using MUSCLE 3

Yvette Pijnenburg

July 27, 2020

Supervisor(s): Pavel Zun (Erasmus MC; IvI, UvA);

Inf

orma

tica

—

Universiteit

v

an

Ams

terd

am

(2)

(3)

Abstract

As cardiovascular diseases are one of the leading causes of death around the world, there are many studies about what the causes are and how to prevent or treat them. But this is not everything. It can also be useful to predict what the body will do after for example a treatment. This means that the patient needs less monitoring, because we almost certainly know how the patient will recover from the treatment.

One of these tools to predict how the body will react is the ISR3D model. This model models the reaction of the immune system after a stent procedure, used to treat stenosis of the arteries. But before this model can be fully operational, we first need to fine tune it. That will be done by doing uncertainty quantification on the model. Prior to this the model first needs to be able to communicate with MUSCLE 3, so that the uncertainty quantification can be done at all. This means that the current version, which communicates with MUSCLE 2, needs to be rewritten. This thesis will document the process of rewriting the model from Java MUSCLE 2 code to C++ MUSCLE 3 code, by documenting the translating process and a few tests to show the results of the translation.

(4)

(5)

2.3 MUSCLE . . . 11 2.3.1 MUSCLE 2 . . . 12 2.3.2 MUSCLE 3 . . . 13 3 My work 15 3.1 Goal . . . 15 3.2 Setup . . . 15 3.3 Implementation . . . 15 3.3.1 Mappers . . . 18 3.3.2 Classes . . . 19 3.4 Testing . . . 20 3.5 Results . . . 21 4 Conclusions 23 Bibliography 25 5 Appendices 27 5.1 Appendix A . . . 27 5.2 Appendix B . . . 31 5.3 Appendix C . . . 36

(6)

(7)

CHAPTER 1

Introduction

Cardiovascular diseases are one of the leading causes of death around the world except Africa (Mendis et al.,2011). One type of cardiovascular disease is coronary artery disease (CAD) and specifically the stenosis, abnormal narrowing, of coronary arteries. This leads to a restricted blood flow to the heart, possibly resulting in a heart attack. This narrowing can be corrected by balloon angioplasty. During this procedure a stent, a mesh-like support, can be placed inside the vessel to keep it open after the vessel lumen is enlarged by the balloon. This procedure also damages the vessel wall, both by stretching it and placing the stent, causing a healing response. This can lead to stent thrombosis, the formation of a blood clot that can block the artery, or in-stent restenosis (ISR), a repeated abnormal narrowing of the artery, due to excessive growth (Jukema et al.,2012).

The probable causes of ISR, as well as possible solutions, were studied in many clinical trials and animal experiments, and also there are some analytical and computational models of it (P. S. Zun et al.,2017).

One of these models is the ISR3D model (P. Zun et al., 2019). This is a multi-scale model of tissue growth and remodelling in arteries after a surgical intervention. ISR3D couples a lattice-based model of flow to an agent-lattice-based model of arterial wall as seen in Figure 2.3. For enabling easy uncertainty quantification and sensitivity analysis of the model, it has to be ported to use MUSCLE 3 (Multiscale Coupling Library and Environment 3) for coupling. The coupling and mapping values from agents to lattice cells and vice versa is done via special modules called mappers.

There are already mappers present in ISR3D, but these are written in Java and MUSCLE 3 does not support Java, but it does support C++. So if the mappers are to be translated to C++ from Java, the model can be coupled with MUSCLE 3 and the uncertainty quantification and sensitivity analysis of the model can be done. Hereby fine tuning the model so the healing process of the stent procedure can be modelled and the doctors know immediately after the procedure if it was a success or that the patient needs to be monitored.

This thesis will document the process of the translation of the mappers from Java MUSCLE 2 code to C++ MUSCLE 3 code, by documenting the translating process and a few tests to show the results of the translation.

(8)

(9)

CHAPTER 2

Theoretical background

2.1 The stenting procedure

Figure 2.1: Balloon

angio-plasty and placement of the

stent. Source: Wikipedia, https:

//commons.wikimedia.org/wiki/File: PTCA stent NIH.gif.

The heart is a muscular organ which pumps

blood through the body. It pumps about 6

litres of blood per minute. The heart

it-self also needs blood to do its work. It gets

this through the coronary arteries. When

suf-fering from coronary artery disease, the

coro-nary arteries are narrowed or blocked, due to a

build-up of plaque, which mostly consists of

fi-brous tissue and lipid cells that infiltrate from the blood stream into the wall (Tarbell et al.,

2014). This can result in an inadequate blood

supply, leaving the heart to become starved of

oxygen and other vital nutrients it needs to

work properly. Possibly resulting in a heart

at-tack.

One of the treatments is balloon angioplasty (Castaneda-Zuniga et al., 1980). This is a minimally invasive pro-cedure to open coronary arteries, which involves a in-flatable balloon-tipped catheter. The catheter is guided through the artery to the blockage. The balloon is po-sitioned in the stenosed segment and then inflated to open up the artery. Usually a stent is placed to keep the artery open after the balloon is deflated and with-drawn.

When placing the stent, the struts of the stent are forced into contact with the vessel wall, hereby widening the vessel diameter as seen in Figure 2.1. Immediately after the procedure, most of the endothe-lium, the inner lining of the blood vessel, is damaged and dysfunctional. The smooth muscle cells (SMCs), responsible for among other things blood flow distribution (Bowens and Parmacek,

2012), change their phenotype, behaviour, from contractile to synthetic and start proliferating (P. S. Zun et al., 2017). Normally the SMCs multiply enough to cover the stent struts in scar tissue, resulting in a recovered endothelium, which in turn suppresses the SMC growth, causing them to switch back to their normal behaviour. However the body’s healing process is not al-ways perfect. Sometimes the growth is too intensive and causes a repeat narrowing of the artery, in-stent restenosis (ISR), see Figure 2.2.

(10)

Figure 2.2: The PCI procedure resulting in ISR. Source: Theranostics,https://www.thno.org/ v04p0175.htm.

2.2 The ISR3D model

The ISR3D model is a three-dimensional multi-scale model of in-stent restenosis in which SMC proliferation and re-endothelialization are fully coupled to blood flow simulations. Figure 2.3 shows the relationships between the main submodels. The ISR3D model is based upon its 2D version, ISR2D (Caiazzo et al.,2011).

ISR is a complex process influenced by factors acting on different scales. A single-scale model would require too much computation time, so a multi-scale modelling approach is the solution. In the ISR3D model the blood flow, the scar tissue and the different layers of the vessels, consisting of the internal elastic lamina (IEL), the media, and the external elastic lamina (EEL), are explicitly modelled. The endothelium is modelled implicitly and the adventitia, the outer layer of tissue surrounding the artery, is assumed to hold the EEL in set position after the stent is deployed (P. S. Zun et al.,2017).

As seen in Figure 2.3, ISR3D is composed of several single-scale submodels, each with their own scale, that are used to model the blood flow inside the stented segment of the artery, the stent deployment and the cell growth, and communicating relevant values to each other.

Figure 2.3: An overview of the multi-scale in-stent restenosis model, illustrating data flow between single-scale components. After simulating stent deployment, post-stenting geometry is used as a starting point for the off-lattice agent-based simulation of neointima, scar tissue, formation. From agent sizes and positions, a voxel-based geometry of the vessel is constructed, which is used to calculate flow in the vessel. The flow produces shear stress on the vessel wall, which is mapped back to the agent-based model to inform the next iteration of neointima formation simulation (P. Zun et al.,2019).

The latest version of ISR3D has been proven to closely match in vivo data (P. Zun et al.,2019). A detailed description of the first ISR3D model can be found in P. S. Zun et al.,2017. The latest version is extended to include extracellular matrix (ECM) production and a modified mechanical

(11)

model for IEL and EEL (P. Zun et al.,2019).

To further improve this model, the goal is to do Uncertainty Quantification (UQ) on it. This has been done on the 2D version of this model (Nikishova et al.,2018). It showed how a variation of parameters affected the simulation results and how important the exact value of a parameter is for the precision of them. As the ISR2D model is very similar, but more lightweight, seeing that it only had five subprocesses in 2D, it is not that simple to do UQ for ISR3D. For a compari-son between the two models see Figure 2.4, where (a) shows the simulated artery and the yet to be placed stent, (b) the stent deployment and (c) the post-deployment SMC proliferation, and Figure 2.5, where (A) shows the vessel directly after stenting and (B) 35 days after stenting.

Figure 2.4: ISR2D initial state, post-deployment state, and re-endothelialization pattern. Orange circles are individual SMCs, red circles are IEL agents, and grey circles make up the stent struts. Green cells are synthetic and proliferating, blue are nitric oxide (NO) inhibited. Note the rupture of the IEL in the vicinity of the strut, and subsequent SMC proliferation in this area (Nikishova et al.,2018).

Figure 2.5: Real-sized stent simulation; (A) right after stenting; (B) 35 days after stenting. The growth happens outside the stented area because we assume that the endothelium is denudated in the entire simulated segment (P. S. Zun et al.,2017).

2.3 MUSCLE

MUSCLE stands for Multi-scale Coupling Library and Environment. In this section we will discuss the two latest versions and their differences in relation to this work.

(12)

2.3.1 MUSCLE 2

The second version of MUSCLE was published in 2014 (Borgdorff et al.,2014) and is the successor of the first version of 2008 (Hegewald et al.,2008), then not yet called MUSCLE. MUSCLE 2 treats single scale models as separate components and makes it possible to couple them, whether they run in one location or multiple. It takes advantage of the fact that submodels can be seen as components in a component-based simulation. Each component tracks their own progress and synchronises with others when they exchange messages.

One of the requirements of MUSCLE 2 is that the submodels are clearly separated, so that a submodel does not dictate how it should be coupled to another submodel. Instead, a submodel sends and receives messages through specified ports, which are then later coupled to each other.

MUSCLE 2 can be divided in three different subdivisions. The API, which submodel code

uses, the coupling scripting environment that tells how the submodels will be coupled, and a runtime environment that runs the multi-scale model on different machines, be they a desktop computer or a supercomputer, or run it on one of these machines. These three subdivisions are all independent from each other, as visualised in Figure 2.6. This means that with one implementation of a submodel, it can be coupled multiple ways before running it.

Figure 2.6: An overview of the communication of MUSCLE 2 and a multi-scale model. Source: MUSCLE 2 Wiki,http://apps.man.poznan.pl/trac/muscle/wiki.

The MUSCLE 2 library consists of Java, C, C++, Python, Fortran APIs and coupling definitions in Ruby. The Java API allows implementing formal Multiscale Modelling Language (MML) constructs like filters and mappers. Because MUSCLE 2 supports multiple languages in a multi-scale model, MML constructs can be used in models that don’t use Java.

The whole documentation of MUSCLE 2 can be found in the paper of Borgdorff et al.,2014and the software can be found athttp://apps.man.poznan.pl/trac/muscle/wiki.

(13)

2.3.2 MUSCLE 3

MUSCLE 3 is the successor to MUSCLE 2. It does the same as MUSCLE 2, it makes creating multi-scale simulations easy, but it also enables efficient UQ of those models by enabling large-scale Monte Carlo experiments as well as advanced semi-intrusive algorithms. It is still being developed by the e-MUSC project of the University of Amsterdam and the Netherlands eScience Center. This means that not all the features are yet implemented and that it is a work in progress. The documentation can be found at https://muscle3.readthedocs.io/en/latest/index.html and the code at https://github.com/multiscale/muscle3(Veen,2020).

MUSCLE 3 uses the Multiscale Modelling and Simulation Language (MMSL) to describe the structure of a multi-scale model. MMSL lets us describe which compute elements, submodels or UQ components for instance, a model consists of, how many instances of each we need, and how they are put together.

Currently MUSCLE 3 supports Python 3, C++ and Fortran, and that is the main technical difficulty for porting models between the two versions for this project. ISR3D has mappers written in Java so that it can be coupled using MUSCLE 2. If we want to do UQ, we need to port ISR3D to MUSCLE 3.

(14)

(15)

CHAPTER 3

My work

3.1 Goal

As stated before in chapter 1, the goal of this project is the documentation of the translation of the mappers of ISR3D from Java MUSCLE 2 code to C++ MUSCLE 3 code. In this chapter we will talk about the setup of the project and the necessary libraries, the implementation of the mappers and the results of the new implementation.

3.2 Setup

As described in chapter 2 we are depending on three different programs to run and test this project. These programs needed to be installed first. Some are easy to install, others need spe-cific versions of their dependencies.

MUSCLE 2 is easy to install. After downloading the directory form the website (https:// github.com/psnc-apps/muscle2) and installing the dependencies with sudo apt-get install build-essential cmake ruby default-jdk, the only command needed to install MUSCLE 2 is sudo ./build.sh in the build folder of the directory.

The source code of ISR3D is located athttps://gitlab.computationalscience.nl/pavel/ISR3Dand for this project the muscle3-port-branch is used. As ISR3D is not widely used, there did not exist an up-to-date installation guide. Some dependencies required specific or older versions than the standard installation version. At the end of the project a new installation guide was created with these cases documented.

MUSCLE 3 is also easy to install and for this project the MUSCLE 3 version used is the version of 10 July 2020 on the develop-branch (https://github.com/multiscale/muscle3/tree/develop).

3.3 Implementation

As described in section 2.3.2 we want to rewrite or translate certain parts of the ISR3D model so that it can communicate with MUSCLE 3 instead of MUSCLE 2. This amounts to the transla-tion of three mapper files from Java MUSCLE 2 code to C++ MUSCLE 3 code and the necessary ISR3D dependency files from Java code to C++ code. As the syntax of Java is similar to the syntax of C/C++, it is not hard to transfer the functionality of Java code to C++ code, for example for-loops and if-statements. A difference is that in Java everything is in a class and the filename is the same as the class in it, while in C++ we do not need to a class if it is not wanted. See Listings 3.1 and 3.2 for some example code with the same result, writing Hello World in a terminal, in Java and C++ respectively.

(16)

1 c l a s s H e l l o W o r l d {

2 p u b l i c s t a t i c v o i d main ( S t r i n g [ ] a r g s ) { 3 System . o ut . p r i n t l n (” H e l l o World ”) ; 4 }

5 }

Listing 3.1: Java code

1 #i n c l u d e 2 3 i n t main ( ) { 4 s t d : : c o u t << ” H e l l o World ” << s t d : : e n d l ; 5 } Listing 3.2: C++ code

On the other hand, even though MUSCLE 2 and 3 do the same things, the API of MUSCLE 2 is very different from MUSCLE 3. For example: to create communication channels as seen in Figure 3.1 in MUSCLE 2, we need to create a variable object ConduitEntrance or Condui-tExit, depending on if the channel is sending or receiving respectively, and can then use them to send or receive data. In MUSCLE 3 we need to create an Instance object with operators, O_F, F_INIT, S, O_I, describing whether the channel is sending or receiving and whether it is at the start or end of the communication or every time step. See Listings 3.3 and 3.4 for a simplified example. 1 p u b l i c c l a s s C e l l M a p p e r f a n I n { 2 // t o SMC 3 p r i v a t e ConduitEntrance<d o u b l e[ ] > M a x S t r e s s S e n d e r ; 4 // from V o x e l i z e r 5 p r i v a t e C o n d u i t E x i t <d o u b l e[ ] > c e l l s E x i t ; 6 p r i v a t e C o n d u i t E x i t <d o u b l e[ ] > domainBoxExit ; 7 p r i v a t e C o n d u i t E x i t c e l l P o s i t i o n E x i t ; 8 p r i v a t e C o n d u i t E x i t domainExit ; 9 p r i v a t e C o n d u i t E x i t s u r f a c e P o s i t i o n E x i t ; 10 11 p r o t e c t e d v o i d a d d P o r t a l s ( ) { 12 // To SMC 13 M a x S t r e s s S e n d e r = addEntrance (” C e l l M a x S t r e s s E n t r a n c e ”, 14 d o u b l e[ ] . c l a s s) ; 15 16 // from V o x e l i z e r 17 c e l l s E x i t = a d d E x i t (” C e l l s E x i t ”, d o u b l e[ ] . c l a s s) ; 18 domainBoxExit = a d d E x i t (” domainBoxExit ”, d o u b l e[ ] . c l a s s) ; 19 c e l l P o s i t i o n E x i t = a d d E x i t (” C e l l P o s i t i o n s E x i t ”, 20 I n t A r r a y L i s t [ ] [ ] [ ] . c l a s s) ; 21 domainExit = a d d E x i t (” DomainExit ”, b y t e[ ] [ ] [ ] . c l a s s) ; 22 s u r f a c e P o s i t i o n E x i t = a d d E x i t (” S u r f a c e P o s i t i o n s E x i t ”, 23 i n t[ ] [ ] . c l a s s) ; 24 } 25 }

(17)

1 /∗ C o l l e c t o r c l a s s ∗/ 2 c l a s s C o l l e c t o r { 3 p r i v a t e: 4 I n s t a n c e& i n s t a n c e ; 5 6 p u b l i c: 7 C o l l e c t o r ( I n s t a n c e& i n s t a n c e ) : 8 i n s t a n c e ( i n s t a n c e ) 9 } ; 10 11 /∗ C r e a t i n g an i n s t a n c e and t h e c o l l e c t o r ∗/ 12 v o i d r u n C o l l e c t o r (i n t a r g c , c h a r∗ a r g v [ ] ) { 13 I n s t a n c e i n s t a n c e ( a r g c , argv , { 14 { Op er at or : : F INIT , {” c e l l s I n ”, 15 ” domainBoxIn ”, 16 ” domainIn ”, 17 ” s u r f a c e P o s i t i o n s I n ”, 18 ” s h e a r S t r e s s I n ”, 19 ” g e o m e t r y I n ”} } , 20 { Op er at or : : O F , {” c e l l M a x S t r e s s O u t ”} } } ; 21 22 C o l l e c t o r c o l l e c t o r ( i n s t a n c e ) ; 23 }

Listing 3.4: C++ MUSCLE 3 code

Figure 3.1: A detailed description of communications, input/output, from the mappers in MUS-CLE 2.

This means that the mappers have to be rewritten and changed instead of simply one-on-one translation. Furthermore, because MUSCLE 3 is still under development, some functionality may not yet work and there is not much reference or examples yet available.

For an overview of all communications between the mappers and submodules in MUSCLE 2 see Figure 3.1. In MUSCLE 3 the submodules blob and drug diffusion are removed, because of deprecation, and with that the thrombus mapper, as seen in Figure 3.2. Despite the removal of those modules, most communications between modules stays the same, save the communications relevant to those modules. This means that only the collector, voxelizer and distributor mapper files need to be translated.

(18)

Figure 3.2: A detailed description of communications, input/output, from the mappers in MUS-CLE 3.

In addition to these three mapper files that need to be translated, there are some other files or functions that need translation from Java to C++ also. Sometimes this is a whole file that defines a variable object, for example the SphereVoxelizer class. Other times it means that there are a few extra lines added in the mapper file as public variables or functions, for example the two helper functions now in the distributor, but original from the Domain class.

Now that we know what needs to be translated, and a little bit how, we will go deeper into these files.

3.3.1 Mappers

CellMapperFanIn/Collector

CellMapperFanIn is as seen in Figure 3.1 the collector. It collects data from the blood flow model, the distributor, or CellMapperFanOut, and the voxelizer and sends data to the SMC module. In the MUSCLE 2 version this file is called CellMapperFanIn, in MUSCLE 3 it is called collector.

This file contains the Collector class and a runCollector function in addition to a main function. The main function calls the runCollector function which in turn creates an instance object like the one in Listing 3.4, calls the Collector class and executes it, by calling the execution function of the Collector class.

The Collector class consists of a constructor, a resetting function and a execution function. The constructor sets the instance it gets form the runCollector function, creates a log file and sets some starting values of variables. The execution function first receives all data from the blood flow model, the distributor and the voxelizer, then unpacks and maps the received values from grid points to SMC agents, logs its progress and checks to see if every grid and data mask have the same size, before sending the resulting data to the SMC and calling the resetting function. The resetting function cleans all variables so that there is no resulting data that can interfere with the next pass.

(19)

CellMapperFanOut/Distributor

CellMapperFanOut is the distributor and collects data from the voxelizer and sends data to the blood flow model and the collector, or CellMapperFanIn. This file also has a different name in the MUSCLE 2 version than in the MUSCLE 3 version, CellMapperFanOut and distributor respectively.

This file contains the Distributor class and a runDistributor function in addition to a main function. The main function calls the runDistributor function which in turn creates an instance object like the one in Listing 3.4, calls the distributor class and executes it, by calling the execute function of the Distributor class.

The Distributor class consists of a constructor, two helper functions, a sending function and an execution function. The constructor sets the instance it gets from the runDistributor function, creates a log file and sets some starting values of variables. The two helper functions determine if the node is a surface node, and checks if the node has an indirect neighbour with a given value, which is used in the surface-finding function. The sending function computes incremental boundaries and sends the result to the blood flow model and the collector. The executing function receives the data from the voxelizer, unpacks and computes it while logging its progress to the log file. When done it clears its data from this run.

Voxelizer

Voxelizer receives data from SMC and sends data to the collector, or CellMapperFanIn, and the distributor, or CellMapperFanOut. It calculates a voxel-based solid domain mask from the coordinates and sizes of the spheres received from the SMC module, before sending the results to the correct receivers. Voxelizer actually consists of two different files. One file calls the Voxelizing class and executes it, while the other contains the class itself. The first file also creates the instance as in Listing 3.4 before it is passed to the class itself.

The Voxelizing class consists of a constructor and an execution function. The constructor sets the instance it gets form the voxelizer function from the first file and creates a log file. The execution function sets the starting values of some variables and sets others by computing the received data from the SMC module before sending the correct results to the correct locations. It also notes major steps and values in its log file.

3.3.2 Classes

Sphere

Sphere describes the logic of a geometric sphere. This file contains a constructor, a operator overloading, getters, setters and some helper functions. The constructor creates a sphere on the basis of three coordinates, a radius and a typeID, which tells us whether the Sphere is a part of the stent or the vessel wall, and in the latter case, its tissue type. The operator that is being overloaded is the =-operator. It now creates a correct sphere object when called. The getters, setters and helper functions are there so that we can access some private variables outside the class, calculate positions of the spheres and check if two spheres are equal to another.

SphereVoxelizer

Where Sphere describes the logic of the sphere, SphereVoxelizer contains the logic for forming a solid voxel domain mask from Sphere-objects, and for mapping values between this voxel domain and the spheres. It contains in total 34 functions, 8 of which are private and the rest are public, and a constructor. The functions are getters, setters or designed to dispose of the values stored in the object. The constructor creates an SphereVoxelizer object from a mesh delta x, a domain, the smoothness radius and a surface detector value to be used during voxel mask calculation. It sets among others these variables and calculates new ones using these variables.

(20)

3.4 Testing

Before we can begin to run simulations for testing the speed of the new version, we first need to test if the mappers are working correctly.

The first step is to see if the translated mappers and classes are building correctly. If they compile we can begin with creating dummy data for testing to see if MUSCLE 3 is working correctly. See Listing 3.5 for how dummy data generating looks in the Voxelizer class. As the other classes, Collector and Distributor, are more primitive than Voxelizer, the test was done in the Voxelizer class. Figure 3.3 shows the generated dummy spheres from Listing 3.5 representation in the 3D voxel domain. 1 // A s p h e r e o b j e c t l o o k s a s f o l l o w s : 2 // [ ( a l i s t o f x , y , z−v a l u e s ) , a r a d i u s , a t y p e I d ] 3 4 // c r e a t e two o b j e c t s 5 a u t o dummyData = Data : : n i l s ( 2 0 ) ; 6 f o r (i n t i = 0 ; i < 2 0 ; i ++) { 7 // s e t t y p e I d 8 i n t 6 4 t t y p e I d = 2 ; 9 // c r e a t e s p h e r e o b j e c t w i t h r e l a t i v e p o s i t i o n s

10 dummyData [ i ] = Data : : l i s t ( Data : : l i s t ( i + 0 . , i + 1 . , i + 1 . ) , 3 . ,

11 t y p e I d ) ;

12 } 13

14 // c r e a t e message w i t h dummy timestamp

15 Message dummyMessage ( 1 . 0 , dummyData ) ;

16 // s e t t o r e c e i v e t h e dummy message a s d e f a u l t i f no c o n n e c t i o n 17 // i s found

18 a u t o msg = i n s t a n c e . r e c e i v e (” c e l l P o s i t i o n s I n ”, dummyMessage ) ; 19 // c o n t i n u e h a n d l i n g normal s t e p s o f t h e v o x e l i z e r

Listing 3.5: Dummy data generation

(21)

The next step is to link a model to the mapper, in this case the receiving data should come from the SMC model, and see if we now get the data from the model instead of the dummy data. If this is the case, the connection between the models works and we can begin with running the complete simulation.

3.5 Results

The simulation used for testing is of a tiny curved vessel of 1mm over a time span of one hour up to 96 hours or four days, with a time step of one hour. The simulation was executed on a laptop with an Intel Core i3-4005U CPU with 4 cores and 4GB RAM. The results of the time it took to run the two different versions, MUSCLE 2 and 3, are seen in Figure 3.4. These simulation times are measured with precalculated blood flow values and the actual computing time without the time it takes to initialise the model.

For more details about the vessel and its parameters, see the isr tinier test vessel -files in the cxa-folder of the project. The .ymmsl -file is for MUSCLE 3, the .cxa.rb-file for MUSCLE 2. Both files have the same parameters, but here we also see the difference between the two MUSCLE versions in how the data is read and how the parameters are specified.

The MUSCLE 2 tests are done with the ISR3D version of 7 May 2020 or muscle2-comm-final -tag, the MUSCLE 3 tests with the ISR3D version of 10 July 2020.

Figure 3.4: Results of the simulation as specified in section 3.5

As we can see in Figure 3.4 both versions take almost the same time to run the simulation. The trend in the figure leads us to believe that for longer simulations MUSCLE 2 will be a little faster than its successor MUSCLE 3.

We also took a look at the total time it took to run each simulation, as the values in Figure 3.4 are of the computing time and do not include the initialisation time of each version. For MUSCLE 2 this was around 22 seconds, while for MUSCLE 3 it was around 12 seconds. If these times are added to the computation time, the two version will be even closer together in runtime.

(22)

Figure 3.5: The individual time per module for one hour of simulation in MUSCLE 3

As we can see in Figure 3.2 the whole model is circular, but it does have a starting and finishing point. The SMC module is initialised first and sends data to the voxelizer, where data is send to both the collector and the distributor. Thereafter sends the distributor its data to the collector and the blood flow module. The blood flow module then sends its data to the collector. Now that the collector has got all its necessary data, it then sends it back to the SMC module, where the last calculations are done before the next iteration can begin.

Both the blood flow model as the SMC module are parallelized, where the blood flow model uses the amount of cores specified when running the simulation, in this case two, the SMC module is parallelized with OpenMP it uses all available cores when calculating the interactions.

The fact that the model follows a clear path per iteration means that we can measure the time each module takes to complete its task. These results are seen in Figure 3.5, where we can clearly see that the SMC module is active from the beginning of the simulation until its end. For some amount of this time it is waiting for input from the other modules. This is not visible in the graph, but we can estimate it by subtracting the time it takes the other modules to complete their tasks. As some modules can run simultaneously it is only an estimate, but it means that the SMC module is actually computing data for about 10 to 13 seconds, instead of the full 23 seconds it takes the simulation.

(23)

CHAPTER 4

Conclusions

As stated before in chapter 1 and section 3.1 the goal of this project was to translate the mappers of the ISR3D model from Java MUSCLE 2 to C++ MUSCLE 3 code and document this process. This process has been described in section 3.3, with the confirmation in sections 3.4 and 3.5 that the new version of the model works.

We also can conclude from section 3.5 that the MUSCLE 2 version of the ISR3D model prob-ably will be slightly faster than the MUSCLE 3 version for longer simulations, and for shorter simulations the difference, if there is any, is within 10% of the total execution time. Due to the limitations of the machine available, it was impossible to test different simulation sizes and longer durations of the simulation.

As said before in section 2.2 the goal is to do UQ on the model to fine tune all parameters, determine the most important input parameters of the model, and to estimate the uncertainties in the output, which arise from the uncertainties in these inputs. This will still be the next step in the ISR3D project, after the MUSCLE 3 mappers are thoroughly tested to work as needed. This process will bring the model one step closer to be eligible for clinical uses in the medical world, for example by performing clinical trials in silico. Hopefully, eventually this will mean that patients with a stent need less frequently monitoring after the procedure, at least in this respect, ISR, of the complications that accompany the procedure. Note that the patient still needs monitoring, because nothing is one hundred percent foolproof, and the ISR3D model only models the ISR complication and not, for example, the formation of blood clots.

(24)

(25)

Bibliography

Borgdorff, J., Mamonski, M., Bosak, B., Kurowski, K., Belgacem, M. B., Chopard, B., Groen, D., Coveney, P. V., & Hoekstra, A. G. (2014). Distributed multiscale computing with muscle 2, the multiscale coupling library and environment. Journal of Computational Science, 5 (5), 719–731.

Bowens, N., & Parmacek, M. S. (2012). Development of the smooth muscle cell lineage, In Muscle. Elsevier.

Caiazzo, A., Evans, D., Falcone, J.-L., Hegewald, J., Lorenz, E., Stahl, B., Wang, D., Bernsdorf, J., Chopard, B., Gunn, J., Et al. (2011). A complex automata approach for in-stent restenosis: Two-dimensional multiscale modelling and simulations. Journal of Computa-tional Science, 2 (1), 9–17.

Castaneda-Zuniga, W., Formanek, A., Tadavarthy, M., Vlodaver, Z., Edwards, J., Zollikofer, C., & Amplatz, K. (1980). The mechanism of balloon angioplasty. Radiology, 135 (3), 565– 571.

Hegewald, J., Krafczyk, M., T¨olke, J., Hoekstra, A., & Chopard, B. (2008). An agent-based coupling platform for complex automata, In International conference on computational science. Springer.

Jukema, J. W., Verschuren, J. J., Ahmed, T. A., & Quax, P. H. (2012). Restenosis after pci. part 1: Pathophysiology and risk factors. Nature Reviews Cardiology, 9 (1), 53.

Mendis, S., Puska, P., Norrving, B., Organization, W. H., Et al. (2011). Global atlas on cardio-vascular disease prevention and control. Geneva: World Health Organization.

Nikishova, A., Veen, L., Zun, P., & Hoekstra, A. G. (2018). Uncertainty quantification of a multiscale model for in-stent restenosis. Cardiovascular engineering and technology, 9 (4), 761–774.

Tarbell, J. M., Shi, Z.-D., Dunn, J., & Jo, H. (2014). Fluid mechanics, arterial disease, and gene expression. Annual review of fluid mechanics, 46.

Veen, L. (2020). Muscle 3: The multiscale coupling library and environment.https://doi.org/10. 5281/zenodo.3764713

Zun, P. S., Anikina, T., Svitenkov, A., & Hoekstra, A. G. (2017). A comparison of fully-coupled 3d in-stent restenosis simulations to in-vivo data. Frontiers in physiology, 8, 284. Zun, P., Narracott, A., Chiastra, C., Gunn, J., & Hoekstra, A. (2019). Location-specific

compari-son between a 3d in-stent restenosis model and micro-ct and histology data from porcine in vivo experiments. Cardiovascular engineering and technology, 10 (4), 568–582.

(26)

(27)

CHAPTER 5

Appendices

5.1 Appendix A

Collector.cpp

1 #i n c l u d e <cmath> 2 #i n c l u d e <c s t d l i b > 3 #i n c l u d e 4 #i n c l u d e <ostream> 5 #i n c l u d e <v e c t o r > 6 #i n c l u d e <s t d e x c e p t > 7 8 #i n c l u d e <l i b m u s c l e / l i b m u s c l e . hpp> 9 #i n c l u d e <ymmsl/ymmsl . hpp> 10 11 #i n c l u d e 12 13 #i n c l u d e ” . . / . . / v o x e l i z e r −cpp / s r c / s p h e r e . hpp” 14 #i n c l u d e ” . . / . . / v o x e l i z e r −cpp / s r c / s p h e r e V o x e l i z e r . hpp” 15 #i n c l u d e ” . . / . . / v o x e l i z e r −cpp / s r c / typeNames . hpp” 16 #i n c l u d e ” . . / . . / v o x e l i z e r −cpp / s r c / i n t P o i n t . hpp” 17 #i n c l u d e ” . . / . . / v o x e l i z e r −cpp / s r c / s u r f a c e D e t e c t o r . hpp” 18 19 u s i n g l i b m u s c l e : : Data ; 20 u s i n g l i b m u s c l e : : DataConstRef ; 21 u s i n g l i b m u s c l e : : I n s t a n c e ; 22 u s i n g l i b m u s c l e : : Message ; 23 u s i n g ymmsl : : O p e r a t o r ; 24 25 u s i n g u t i l : : l o g g e r ; 26 27 c l a s s C o l l e c t o r { 28 p r i v a t e: 29 d o u b l e d e l t a X ; 30 d o u b l e t c u r ; 31 I n s t a n c e& i n s t a n c e ; 32 c o n s t s t d : : s t r i n g outPath ; 33 S p h e r e V o x e l i z e r ∗ s p h e r e v o x = n u l l p t r ; 34 35 s t d : : v e c t o r <d o u bl e> p r e v i o u s S h e a r P a l a b o s D a t a ;

(28)

36 s t d : : v e c t o r p r e v i o u s S o l i d ; 37 s t d : : v e c t o r <d o u bl e> drug ;

38 s t d : : v e c t o r <d o u bl e> n e o i n t i m a ;

39 s t d : : v e c t o r <d o u bl e> c e l l S h e a r p a l a b o s ; 40 // {minX , minY , minZ , maxX , maxY , maxZ} 41 s t d : : v e c t o r <d o u bl e> domainBoundingBox ; 42 s t d : : v e c t o r <s i z e t > b o x S i z e ; // number o f v o x e l s a l o n g X, Y, Z 43 44 p u b l i c: 45 v o i d r e s e t ( ) ; 46 v o i d e x e c u t e ( ) ; 47 48 C o l l e c t o r ( I n s t a n c e& i n s t a n c e ) : 49 i n s t a n c e ( i n s t a n c e ) , 50 outPath ( i n s t a n c e . g e t s e t t i n g a s <s t d : : s t r i n g >( 51 ” m u s c l e d a t a o u t d i r ”) ) { 52 l o g g e r ( ) . s e t F i l e ( outPath + ” c o l l e c t o r . l o g ”) ; 53 d e l t a X = i n s t a n c e . g e t s e t t i n g a s <d o u bl e>(” f l o w d x ”) ; 54 l o g g e r ( ) << ” C o l l e c t o r i n i t i a l i s e d ” << s t d : : e n d l ; 55 } 56 } ; 57 58 v o i d C o l l e c t o r : : e x e c u t e ( ) { 59 a u t o d msg = i n s t a n c e . r e c e i v e (” domainBoxIn ”) ; 60 i f ( s p h e r e v o x == n u l l p t r ) { 61 a u t o domainPtr = d msg . d a t a ( ) . e l e m e n t s <d o u bl e> ( ) ;

62 domainBoundingBox . a s s i g n ( domainPtr , domainPtr +

63 d msg . d a t a ( ) . s i z e ( ) ) ;

64

65 l o g g e r ( ) << ” C r e a t i n g v o x e l i z e r ” << s t d : : e n d l ; 66 N o S u r f a c e s u r f a c e D e t e c t o r ;

67 s p h e r e v o x = new S p h e r e V o x e l i z e r ( deltaX , domainBoundingBox ,

68 0 . , s u r f a c e D e t e c t o r ) ; 69 } 70 71 // t o r e c e i v e c e l l l i s t s from SMC 72 l o g g e r ( ) << ” R e c e i v i n g c e l l p o s i t i o n ” << s t d : : e n d l ; 73 a u t o c msg = i n s t a n c e . r e c e i v e (” c e l l s I n ”) ; 74 75 s t d : : v e c t o r <Sphere> s p h e r e s ; 76 a u t o s p h e r e s D a t a = c msg . d a t a ( ) ; 77 // . . . and p a r s e a l l r e c e i v e d d a t a f i e l d s 78 f o r ( s i z e t i = 0 ; i < s p h e r e s D a t a . s i z e ( ) ; ++i ) { 79 d o u b l e x = s p h e r e s D a t a [ i ] [ 0 ] [ 0 ] . as<d o u bl e> ( ) ; 80 d o u b l e y = s p h e r e s D a t a [ i ] [ 0 ] [ 1 ] . as<d o u bl e> ( ) ; 81 d o u b l e z = s p h e r e s D a t a [ i ] [ 0 ] [ 2 ] . as<d o u bl e> ( ) ; 82 d o u b l e r = s p h e r e s D a t a [ i ] [ 1 ] . as<d o u bl e> ( ) ; 83 // t y p e o f t i s s u e ; s o f a r , i t ’ s smc o r s t e n t 84 i n t typeID = s p h e r e s D a t a [ i ] [ 2 ] . as ( ) ; 85 86 // A new S p h e r e o b j e c t i s made f o r e a c h r e c e i v e d i t e m 87 s p h e r e s . p u s h b a c k ( S p h e r e ( x , y , z , r , typeID ) ) ; 88 } 89 l o g g e r ( ) << s p h e r e s . s i z e ( ) << ” c e l l s r e c e i v e d . ” << s t d : : e n d l ; 90 s p h e r e v o x −>s e t S h a p e s ( s p h e r e s ) ; 91

(29)

92 // r e c e i v e mesh domain 93 // s i m i l a r t o d i s t r i b u t o r , can be moved t o a f u n c t i o n 94 a u t o domainObs = i n s t a n c e . r e c e i v e (” domainIn ”) ; 95 t c u r = domainObs . timestamp ( ) ; 96 l o g g e r ( ) << ”Domain r e c e i v e d f o r t i m e ” << t c u r << s t d : : e n d l ; 97 98 s i z e t d o S i z e = domainObs . d a t a ( ) . s i z e ( ) ; 99 b o x S i z e = domainObs . d a t a ( ) . s h a p e ( ) ; 100 l o g g e r ( ) << ”Domain s i z e ” << d o S i z e << ” e l e m e n t s , ” << 101 b o x S i z e [ 0 ] << ” by ” < aggregateDomain ;

105 a u t o domainObsPtr = domainObs . d a t a ( ) . e l e m e n t s ( ) ;

106 aggregateDomain . a s s i g n ( domainObsPtr , domainObsPtr + d o S i z e ) ;

107 s p h e r e v o x −>setMaskedDomain ( aggregateDomain ) ; 108 l o g g e r ( ) << ” S p h e r e v o x e l i z e r domain a s s i g n e d ” << s t d : : e n d l ; 109 110 // can ’ t r e c e i v e r e f e r e n c e s 111 // a u t o cp msg = i n s t a n c e . r e c e i v e ( ” c e l l P o s i t i o n E x i t ” ) ; 112 // s p h e r e v o x −>s e t S h a p e R e f e r e n c e ( cp . d a t a ( ) ) ; 113 // r e b u i l d v o x e l i z a t i o n l o c a l l y i n s t e a d 114 s p h e r e v o x −>v o x e l i z e N o S u r f a c e ( ) ; 115 116 a u t o sp msg = i n s t a n c e . r e c e i v e (” s u r f a c e P o s i t i o n s I n ”) ; 117 a u t o spData = sp msg . d a t a ( ) ; 118 s t d : : v e c t o r s u r f a c e P o i n t s ; 119 f o r ( s i z e t i = 0 ; i < spData . s i z e ( ) ; ++i ) { 120 s u r f a c e P o i n t s . p u s h b a c k ( i n t p o i n t ( spData [ i ] [ 0 ] . as ( ) , 121 spData [ i ] [ 1 ] . as ( ) , 122 spData [ i ] [ 2 ] . as ( ) ) 123 ) ; 124 } 125 126 l o g g e r ( ) << ” R e c e i v i n g BF d a t a ” << s t d : : e n d l ; 127 s p h e r e v o x −>s e t S u r f a c e ( s u r f a c e P o i n t s ) ; 128 129 // r e c e i v e s u r f a c e v a l u e s from BF 130 l o g g e r ( ) << ” R e c e i v i n g BF d a t a ” << s t d : : e n d l ; 131 132 a u t o c o o r d m s g = i n s t a n c e . r e c e i v e (” g e o m e t r y I n ”) ; 133 s t d : : v e c t o r i s S o l i d ; 134 a u t o s o l P t r = c o o r d m s g . d a t a ( ) . e l e m e n t s ( ) ; 135 i s S o l i d . a s s i g n ( s o l P t r , s o l P t r + c o o r d m s g . d a t a ( ) . s i z e ( ) ) ; 136 137 a u t o p a l a b o s m s g = i n s t a n c e . r e c e i v e (” s h e a r S t r e s s I n ”) ; 138 s t d : : v e c t o r <d o u bl e> s h e a r p a l a b o s D a t a ; 139 a u t o s h e a r P t r = p a l a b o s m s g . d a t a ( ) . e l e m e n t s <d o u bl e> ( ) ; 140 s h e a r p a l a b o s D a t a . a s s i g n ( s h e a r P t r , s h e a r P t r + 141 p a l a b o s m s g . d a t a ( ) . s i z e ( ) ) ; 142 143 i f ( i s S o l i d . s i z e ( ) != 0 ) { 144 l o g g e r ( ) << ”New BF r e s u l t s ” << s t d : : e n d l ; 145 // S t o r e t h e changed v a l u e s 146 p r e v i o u s S o l i d = i s S o l i d ; 147 p r e v i o u s S h e a r P a l a b o s D a t a = s h e a r p a l a b o s D a t a ;

(30)

148 } e l s e { 149 // I f no change happened , u s e p r e v i o u s v a l u e s 150 l o g g e r ( ) << ”No change i n BF r e s u l t s , u s i n g p r e v i o u s r e s u l t s ” 151 << s t d : : e n d l ; 152 i s S o l i d = p r e v i o u s S o l i d ; 153 s h e a r p a l a b o s D a t a = p r e v i o u s S h e a r P a l a b o s D a t a ; 154 } 155 156 i f ( i s S o l i d . s i z e ( ) != s h e a r p a l a b o s D a t a . s i z e ( ) ) { 157 throw s t d : : r u n t i m e e r r o r ( s t d : : s t r i n g (” S i z e o f d a t a mask ( ”) + 158 s t d : : t o s t r i n g ( i s S o l i d . s i z e ( ) ) + 159 ” ) i s n o t e q u a l t o t h e s i z e o f t h e d a t a ( ” + 160 s t d : : t o s t r i n g ( s h e a r p a l a b o s D a t a . s i z e ( ) ) + ” ) ”) ; 161 } 162 i f ( i s S o l i d . s i z e ( ) != b o x S i z e [X] ∗ b o x S i z e [Y] ∗ b o x S i z e [ Z ] ) { 163 throw s t d : : r u n t i m e e r r o r ( s t d : : s t r i n g (” S i z e o f d a t a mask ( ”) + 164 s t d : : t o s t r i n g ( i s S o l i d . s i z e ( ) ) + 165 ” ) i s n o t e q u a l t o t h e s i z e o f t h e g r i d ( ” + 166 s t d : : t o s t r i n g ( ( b o x S i z e [X] ∗ b o x S i z e [Y] ∗ b o x S i z e [ Z ] ) ) + ” ) ”) ; 167 } 168 169 l o g g e r ( ) << ” Mapping w a l l s h e a r s t r e s s ” << s t d : : e n d l ; 170 // / c o n v e r t t o i n t 3 2 t f o r v o x e l i z e r 171 s t d : : v e c t o r i s S o l i d F o r V o x ( i s S o l i d . b e g i n ( ) , 172 i s S o l i d . end ( ) ) ; 173 c e l l S h e a r p a l a b o s = s p h e r e v o x −>mapSurfaceAverage ( s h e a r p a l a b o s D a t a , 174 i s S o l i d F o r V o x ) ; 175 176 l o g g e r ( ) << ”WSS s i z e : ” << c e l l S h e a r p a l a b o s . s i z e ( ) << s t d : : e n d l ; 177 178 // s e n d mapped d a t a t o SMC 179 l o g g e r ( ) << ” S e n d i n g d a t a t o SMC. ” << s t d : : e n d l ; 180 a u t o c e l l r e s u l t = Data : : g r i d ( c e l l S h e a r p a l a b o s . d a t a ( ) , 181 { c e l l S h e a r p a l a b o s . s i z e ( ) } ) ; 182 i n s t a n c e . s e n d (” c e l l M a x S t r e s s O u t ”, Message ( t c u r , c e l l r e s u l t ) ) ; 183 l o g g e r ( ) << ” Data t o SMC s e n t . ” << s t d : : e n d l ; 184 185 r e s e t ( ) ; 186 187 } 188 189 v o i d C o l l e c t o r : : r e s e t ( ) { 190 t c u r = 0 . ; 191 c e l l S h e a r p a l a b o s . c l e a r ( ) ; 192 s p h e r e v o x −>d i s p o s e ( ) ; 193 } 194 195 v o i d r u n C o l l e c t o r (i n t a r g c , c h a r ∗ a r g v [ ] ) { 196 I n s t a n c e i n s t a n c e ( a r g c , argv , { 197 { Op er at or : : F INIT , {” c e l l s I n ”, 198 ” domainBoxIn ”, 199 ” domainIn ”, 200 ” s u r f a c e P o s i t i o n s I n ”, 201 ” s h e a r S t r e s s I n ”, 202 ” g e o m e t r y I n ”} } , 203 { Op er at or : : O F , {” c e l l M a x S t r e s s O u t ”} } } ) ;

(31)

204 205 C o l l e c t o r c o l l e c t o r ( i n s t a n c e ) ; 206 207 w h i l e ( i n s t a n c e . r e u s e i n s t a n c e ( ) ) { 208 c o l l e c t o r . e x e c u t e ( ) ; 209 } 210 } 211 212 i n t main (i n t a r g c , c h a r ∗ a r g v [ ] ) { 213 r u n C o l l e c t o r ( a r g c , a r g v ) ; 214 r e t u r n EXIT SUCCESS ; 215 }

5.2 Appendix B

Voxelizer.cpp

1 #i n c l u d e <c s t d l i b > 2 #i n c l u d e < l i s t > 3 #i n c l u d e <l i b m u s c l e / l i b m u s c l e . hpp> 4 #i n c l u d e <ymmsl/ymmsl . hpp> 5 #i n c l u d e ” v o x e l . hpp” 6 7 u s i n g l i b m u s c l e : : Data ; 8 u s i n g l i b m u s c l e : : DataConstRef ; 9 u s i n g l i b m u s c l e : : I n s t a n c e ; 10 u s i n g l i b m u s c l e : : Message ; 11 u s i n g ymmsl : : O p e r a t o r ; 12 13 v o i d v o x e l i z e r (i n t a r g c , c h a r ∗ a r g v [ ] ) { 14 I n s t a n c e i n s t a n c e ( a r g c , argv , { 15 { Op er at or : : F INIT , {” c e l l P o s i t i o n s I n ”} } , 16 { Op er at or : : O F , {” domainOut1 ”, 17 ” domainOut2 ”, 18 ” c e l l s O u t ”, 19 ” s u r f a c e P o s i t i o n O u t ”, 20 ” domainBoxOut ”}} 21 } ) ; 22 23 V o x e l i z i n g M a p p e r v o x e l i z i n g M a p p e r ( i n s t a n c e ) ; 24 25 w h i l e ( i n s t a n c e . r e u s e i n s t a n c e ( ) ) { 26 v o x e l i z i n g M a p p e r . e x e c u t e ( ) ; 27 } 28 } 29 30 i n t main (i n t a r g c , c h a r ∗ a r g v [ ] ) { 31 v o x e l i z e r ( a r g c , a r g v ) ; 32 r e t u r n EXIT SUCCESS ; 33 }

Voxel.hpp

1 #i f n d e f VOXEL HPP 2 #d e f i n e VOXEL HPP 3

(32)

4 #i n c l u d e <cmath> 5 #i n c l u d e <c s t d l i b > 6 #i n c l u d e 7 #i n c l u d e <ostream> 8 #i n c l u d e <v e c t o r > 9 #i n c l u d e <l i m i t s > 10 #i n c l u d e < l i s t > 11 #i n c l u d e <c s t d i o > 12 #i n c l u d e <a r r a y > 13 #i n c l u d e <s t d e x c e p t > 14 15 #i n c l u d e <l i b m u s c l e / l i b m u s c l e . hpp> 16 #i n c l u d e <ymmsl/ymmsl . hpp> 17 18 #i n c l u d e 19 20 #i n c l u d e ” typeNames . hpp” 21 #i n c l u d e ” s p h e r e V o x e l i z e r . hpp” 22 23 u s i n g l i b m u s c l e : : Data ; 24 u s i n g l i b m u s c l e : : DataConstRef ; 25 u s i n g l i b m u s c l e : : I n s t a n c e ; 26 u s i n g l i b m u s c l e : : Message ; 27 u s i n g ymmsl : : O p e r a t o r ; 28 29 u s i n g u t i l : : l o g g e r ; 30 31 c l a s s V o x e l i z i n g M a p p e r { 32 p r i v a t e: 33 d o u b l e l e n g t h ; 34 d o u b l e d e l t a X ; 35 d o u b l e s m o o t h r a d i u s ; 36 i n t s u r f a c e i t e r ; 37 I n s t a n c e& i n s t a n c e ; 38 c o n s t s t d : : s t r i n g outPath ; 39 40 p u b l i c: 41 V o x e l i z i n g M a p p e r ( I n s t a n c e& i n s t a n c e ) : 42 i n s t a n c e ( i n s t a n c e ) , 43 outPath ( i n s t a n c e . g e t s e t t i n g a s <s t d : : s t r i n g >( 44 ” m u s c l e d a t a o u t d i r ”) ) { 45 l o g g e r ( ) . s e t F i l e ( outPath + ” v o x e l . l o g ”) ; 46 } 47 48 v o i d e x e c u t e ( ) ; 49 50 } ; 51 52 v o i d V o x e l i z i n g M a p p e r : : e x e c u t e ( ) { 53 s t a t i c i n t i t e r a t i o n = 0 ; 54 s t a t i c s t d : : v e c t o r <d o u bl e> domainBoundingBox ( 6 ) ; 55 s t a t i c i n t d e t e c t o r Y = 0 ; 56 s t a t i c i n t d e t e c t o r Z = 0 ; 57 58 // A v e c t o r f o r a l l S p h e r e o b j e c t s we r e c e i v e t h i s i t e r a t i o n 59 s t d : : v e c t o r <Sphere> s p h e r e s ;

(33)

60 61 // Read p a r a m e t e r s 62 d e l t a X = i n s t a n c e . g e t s e t t i n g a s <d o u bl e>(” f l o w d x ”) ; 63 s m o o t h r a d i u s = i n s t a n c e . g e t s e t t i n g a s <d o u bl e>(” s m o o t h R a d i u s ”) ; 64 s u r f a c e i t e r = i n s t a n c e . g e t s e t t i n g a s (” s u r f a c e i t e r ”) ; 65 66 l o g g e r ( ) << ” R e c e i v i n g c e l l p o s i t i o n d a t a from AB SMC submodel , 67 i t e r a t i o n ” << i t e r a t i o n << s t d : : e n d l ; 68 // R e c e i v e message . . . 69 a u t o dummyData = Data : : n i l s ( 2 0 ) ; 70 f o r (i n t i = 0 ; i < 2 0 ; i ++) { 71 // p o i n t t pos = { 1 , 2 , 3 } ; 72 i n t 6 4 t t y p e I d = 2 ;

73 dummyData [ i ] = Data : : l i s t ( Data : : l i s t ( i + 0 . , i + 1 . , i + 1 . ) ,

74 3 . , t y p e I d ) ;

75 }

76 Message dummyMessage ( 1 . 0 , dummyData ) ;

77 a u t o msg = i n s t a n c e . r e c e i v e (” c e l l P o s i t i o n s I n ”, dummyMessage ) ; 78 l o g g e r ( ) << ” R e c e i v e d . . . ” << s t d : : e n d l ; 79 // a u t o msg = i n s t a n c e . r e c e i v e ( ” c e l l P o s i t i o n s I n ” ) ; 80 a u t o s p h e r e s D a t a = msg . d a t a ( ) ; 81 // . . . and p a r s e a l l r e c e i v e d d a t a f i e l d s 82 f o r ( s i z e t i = 0 ; i < s p h e r e s D a t a . s i z e ( ) ; ++i ) { 83 d o u b l e x = s p h e r e s D a t a [ i ] [ 0 ] [ 0 ] . as<d o u bl e> ( ) ; 84 d o u b l e y = s p h e r e s D a t a [ i ] [ 0 ] [ 1 ] . as<d o u bl e> ( ) ; 85 d o u b l e z = s p h e r e s D a t a [ i ] [ 0 ] [ 2 ] . as<d o u bl e> ( ) ; 86 d o u b l e r = s p h e r e s D a t a [ i ] [ 1 ] . as<d o u bl e> ( ) ; 87 // t y p e o f t i s s u e ; s o f a r , i t ’ s smc o r s t e n t 88 i n t typeID = s p h e r e s D a t a [ i ] [ 2 ] . as ( ) ; 89 90 // A new S p h e r e o b j e c t i s made f o r e a c h r e c e i v e d i t e m 91 s p h e r e s . p u s h b a c k ( S p h e r e ( x , y , z , r , typeID ) ) ; 92 } 93 94 l o g g e r ( ) << s p h e r e s . s i z e ( ) << ” c e l l s r e c e i v e d . ” << s t d : : e n d l ; 95 96 // make a v o x e l geometry 97 l o g g e r ( ) << ” V o x e l i z i n g ” << s t d : : e n d l ; 98 // assuming t h e o u t e r w a l l i s s t a t i c , we o n l y have t o i n i t i a l i z e 99 // S u r f a c e D e t e c t o r o n c e a t t h e 1 s t s t e p 100 i f ( i t e r a t i o n == 0 ) {

101 d o u b l e minX , minY , minZ ;

102 minX = minY = minZ = s t d : : n u m e r i c l i m i t s <d o u bl e> : : max ( ) ; 103

104 d o u b l e maxX , maxY , maxZ ;

105 maxX = maxY = maxZ = s t d : : n u m e r i c l i m i t s <d o u bl e> : : l o w e s t ( ) ; 106

107 //A bounding box i s c o n s t r u c t e d which i n c l u d e s a l l c e l l s . 108 // I t i s assumed t h a t t h e i n l e t and o u t l e t a r e a l r e a d y 109 // p e r p e n d i c u l a r t o X a x i s and t h e r e i s no l e a k a g e 110 // Look f o r max/min c e l l c o o r d i n a t e s 111 f o r ( S p h e r e& s : s p h e r e s ) { 112 i f ( s .X <= minX ) { 113 minX = s .X; 114 } 115 i f ( s .X >= maxX) {

(34)

116 maxX = s .X; 117 } 118 i f ( s .Y <= minY ) { 119 minY = s .Y; 120 } 121 i f ( s .Y >= maxY) { 122 maxY = s .Y; 123 } 124 i f ( s . Z <= minZ ) { 125 minZ = s . Z ; 126 } 127 i f ( s . Z >= maxZ ) { 128 maxZ = s . Z ; 129 } 130 } 131

132 // add padding around t h e v e s s e l and a d j u s t X−a x i s bounds

133 maxX −= d e l t a X ; 134 minX += d e l t a X ; 135 maxY += 2 ∗ d e l t a X ; 136 minY −= 2 ∗ d e l t a X ; 137 maxZ += 2 ∗ d e l t a X ; 138 minZ −= 2 ∗ d e l t a X ;

139 l o g g e r ( ) << ”Domain : from ( ” << minX << ” , ” << minY << ” , ” << 140 minZ << ” ) t o ( ” << maxX << ” , ” << maxY << ” , ” << maxZ << ” ) . ” 141 << s t d : : e n d l ;

142 domainBoundingBox = {minX , minY , minZ , maxX , maxY , maxZ } ;

143

144 // Look f o r t h e min and max o f t h e l o w e s t −X s l i c e

145 d o u b l e lowestXminY = s t d : : n u m e r i c l i m i t s <d o u bl e> : : max ( ) ; 146 d o u b l e lowestXminZ = s t d : : n u m e r i c l i m i t s <d o u bl e> : : max ( ) ; 147 d o u b l e lowestXmaxY = s t d : : n u m e r i c l i m i t s <d o u bl e> : : l o w e s t ( ) ; 148 d o u b l e lowestXmaxZ = s t d : : n u m e r i c l i m i t s <d o u bl e> : : l o w e s t ( ) ; 149 150 f o r ( S p h e r e& s : s p h e r e s ) { 151 i f ( s .X <= minX + 2 ∗ d e l t a X ) { 152 i f ( s .Y <= lowestXminY ) { 153 lowestXminY = s .Y; 154 } 155 i f ( s .Y >= lowestXmaxY ) { 156 lowestXmaxY = s .Y; 157 } 158 i f ( s . Z <= lowestXminZ ) { 159 lowestXminZ = s . Z ; 160 } 161 i f ( s . Z >= lowestXmaxZ ) { 162 lowestXmaxZ = s . Z ; 163 } 164 } 165 } 166 d e t e c t o r Y = s t a t i c c a s t(( lowestXmaxY+lowestXminY −2∗minY ) 167 / ( 2 . ∗ d e l t a X ) ) ; 168 d e t e c t o r Z = s t a t i c c a s t(( lowestXmaxZ+lowestXminZ −2∗minZ ) 169 / ( 2 . ∗ d e l t a X ) ) ; 170 } 171

(35)

172 // I n i t i a l i z e p e r c o l a t i o n d e t e c t o r w i t h a node t h a t ’ s d e f i n i t e l y 173 // i n s i d e t h e lumen 174 P e r c o l a t i o n D e t e c t o r s u r f a c e D e t e c t o r ( i n t p o i n t ( 0 , d e t e c t o r Y , 175 d e t e c t o r Z ) ) ; 176 177 l o g g e r ( ) << ” C r e a t i n g new sphereVox ” << s t d : : e n d l ; 178 S p h e r e V o x e l i z e r sphereVox ( deltaX , domainBoundingBox ,

179 s m o o t h r a d i u s , s u r f a c e D e t e c t o r ) ; 180 181 l o g g e r ( ) << ” S e t v a l u e s o f sphereVox ” << s t d : : e n d l ; 182 183 sphereVox . s e t S h a p e s ( s p h e r e s ) ; 184 l o g g e r ( ) << ” Shapes s e t ” << s t d : : e n d l ; 185 sphereVox . v o x e l i z e ( ) ; 186 l o g g e r ( ) << ” V o x e l i z e d ” << s t d : : e n d l ; 187 sphereVox . p e r c o l a t e O u t s i d e ( ) ; 188 l o g g e r ( ) << ” P r e l o c a t e d o u t s i d e ” << s t d : : e n d l ; 189 190 // s e t s t e n t l a t t i c e p o i n t s from ” s o l i d ” t o ” s t r u t ” 191 s t d : : v e c t o r aggregateDomain ;

192 aggregateDomain = sphereVox . getSmoothDomain ( ) ;

193 l o g g e r ( ) << ”Smooth domain c a l c u l a t e d ” << s t d : : e n d l ; 194 195 d o u b l e t c u r = msg . timestamp ( ) ; 196 197 a u t o b o x s i z e r e s u l t = Data : : g r i d ( domainBoundingBox . d a t a ( ) , 198 { domainBoundingBox . s i z e ( ) } , {” xyz ”} ) ; 199 i n s t a n c e . s e n d (” domainBoxOut ”, Message ( t c u r , b o x s i z e r e s u l t ) ) ; 200 201 // Count n o d e s by type , i n i t a l l t y p e s w i t h z e r o e s 202 s t d : : a r r a y < s i z e t , 256> numNodes ; 203 numNodes . f i l l ( 0 ) ; 204 f o r (a u t o v a l : aggregateDomain ) { 205 i f ( v a l >= 2 5 6 ) {

206 throw s t d : : r u n t i m e e r r o r (”More than 255 t i s s u e t y p e s a r e

207 n o t s u p p o r t e d ”) ; 208 } 209 numNodes [ v a l ]++; 210 } 211 212 l o g g e r ( ) << ” T o t a l n o d e s : ” << ( aggregateDomain . s i z e ( ) ) << 213 ” , f l u i d n o d e s : ” << numNodes [ FLUID ] << ” , s o l i d : ” <<

214 numNodes [ SOLID ] << ” , s t r u t : ” << numNodes [STRUT] << ” ,

215 s t a t i c n o d e s : ” << numNodes [ STATIC ] << s t d : : e n d l ;

216

217 a u t o a D r e s u l t = Data : : g r i d ( aggregateDomain . d a t a ( ) , { 218 sphereVox . s i z e X , sphereVox . s i z e Y , sphereVox . s i z e Z } ,

219 {”x”, ”y”, ” z ”} ) ; 220 221 i n s t a n c e . s e n d (” domainOut1 ”, Message ( t c u r , a D r e s u l t ) ) ; 222 i n s t a n c e . s e n d (” domainOut2 ”, Message ( t c u r , a D r e s u l t ) ) ; 223 l o g g e r ( ) << ” DomainEntrances s e n t ” << s t d : : e n d l ; 224 225 i f( ( s u r f a c e i t e r != 0 ) && ( i t e r a t i o n % s u r f a c e i t e r == 0 ) ) { 226 c h a r s t r [ 7 ] ; 227 s n p r i n t f ( s t r , 7 , ”%06d”, i t e r a t i o n ) ;

(36)

228 sphereVox . w r i t e S u r f a c e T o F i l e ( outPath + ” / s u r f a c e p o i n t s ” + 229 s t d : : s t r i n g ( s t r ) + ” . xyz ”) ; 230 } 231 i n s t a n c e . s e n d (” c e l l s O u t ”, Message ( t c u r , s p h e r e s D a t a ) ) ; 232 233 s t d : : v e c t o r s u r f = sphereVox . g e t S u r f a c e ( ) ; 234 a u t o s u r f r e s u l t = Data : : n i l s ( s u r f . s i z e ( ) ) ; 235 f o r ( s t d : : s i z e t i = 0 ; i < s u r f . s i z e ( ) ; ++i ) { 236 s u r f r e s u l t [ i ] = Data : : l i s t ( s u r f [ i ] [ 0 ] , s u r f [ i ] [ 1 ] , 237 s u r f [ i ] [ 2 ] ) ; 238 } 239 i n s t a n c e . s e n d (” s u r f a c e P o s i t i o n O u t ”, Message ( t c u r , s u r f r e s u l t ) ) ; 240 l o g g e r ( ) << ” S e n t s u r f a c e p o s i t i o n s t o t h e c o l l e c t o r ” << s t d : : e n d l ; 241 242 i t e r a t i o n ++; 243 l o g g e r ( ) << ”Out o f c e l l mapper l o o p ” << s t d : : e n d l ; 244 } 245 #e n d i f

5.3 Appendix C

Distributor.cpp

1 #i n c l u d e <cmath> 2 #i n c l u d e <c s t d l i b > 3 #i n c l u d e 4 #i n c l u d e <ostream> 5 #i n c l u d e <v e c t o r > 6 7 #i n c l u d e <l i b m u s c l e / l i b m u s c l e . hpp> 8 #i n c l u d e <ymmsl/ymmsl . hpp> 9 10 #i n c l u d e 11 #i n c l u d e ” . . / . . / v o x e l i z e r −cpp / s r c / typeNames . hpp” 12 13 u s i n g l i b m u s c l e : : Data ; 14 u s i n g l i b m u s c l e : : DataConstRef ; 15 u s i n g l i b m u s c l e : : I n s t a n c e ; 16 u s i n g l i b m u s c l e : : Message ; 17 u s i n g ymmsl : : O p e r a t o r ; 18 19 u s i n g u t i l : : l o g g e r ; 20 21 c l a s s D i s t r i b u t o r { 22 p r i v a t e: 23 d o u b l e t c u r ; 24 25 v o i d c o m p u t e A n d S e n d I n c r e m e n t a l B o u n d a r i e s ( 26 c o n s t s t d : : v e c t o r & domain , 27 s t d : : v e c t o r & p r e v i o u s S o l i d ) ; 28 29 I n s t a n c e& i n s t a n c e ; 30 c o n s t s t d : : s t r i n g outPath ; 31 32 s t d : : v e c t o r p r e v i o u s S o l i d ; 33 s t d : : v e c t o r aggregateDomain ;

(37)

34 s t d : : v e c t o r <s i z e t > b o x S i z e ; 35 36 p u b l i c: 37 38 v o i d e x e c u t e ( ) ; 39 40 D i s t r i b u t o r ( I n s t a n c e& i n s t a n c e ) : 41 i n s t a n c e ( i n s t a n c e ) , 42 outPath ( i n s t a n c e . g e t s e t t i n g a s <s t d : : s t r i n g >( 43 ” m u s c l e d a t a o u t d i r ”) ) { 44 l o g g e r ( ) . s e t F i l e ( outPath + ” d i s t r i b u t o r . l o g ”) ; 45 46 p r e v i o u s S o l i d . c l e a r ( ) ; 47 aggregateDomain . c l e a r ( ) ; 48 b o x S i z e = { 0 , 0 , 0 } ; 49 l o g g e r ( ) << ” I n i t i a l i z i n g D i s t r i b u t o r ” << s t d : : e n d l ; 50 } 51 52 b o o l i s S u r f a c e N o d e (i n t i , i n t j , i n t k ,

53 c o n s t s t d : : v e c t o r & nodeTypes ) c o n s t {

54 r e t u r n ( nodeTypes [ i ∗ b o x S i z e [Y] ∗ b o x S i z e [ Z ] + j ∗ b o x S i z e [ Z ] + k ] 55 != FLUID && h a s I n d i r e c t N e i g h b o r W i t h V a l u e ( nodeTypes , i ,

56 j , k , FLUID ) ) ;

57 } 58

59 b o o l h a s I n d i r e c t N e i g h b o r W i t h V a l u e (

60 c o n s t s t d : : v e c t o r & domain , i n t i , i n t j , i n t k ,

61 i n t 3 2 t v a l u e ) c o n s t { 62 f o r (i n t i i = i − 1 ; i i <= i + 1 ; i i ++) { 63 i f ( i i == −1 | | i i == b o x S i z e [X ] ) c o n t i n u e; 64 f o r (i n t j j = j − 1 ; j j <= j + 1 ; j j ++) { 65 i f ( j j == −1 | | j j == b o x S i z e [Y ] ) c o n t i n u e; 66 f o r (i n t kk = k − 1 ; kk <= k + 1 ; kk++) { 67 i f ( kk == −1 | | kk == b o x S i z e [ Z ] | | 68 ( i i == i && j j == j && kk == k ) ) c o n t i n u e; 69 i f ( domain [ i i ∗ b o x S i z e [Y] ∗ b o x S i z e [ Z ] + 70 j j ∗ b o x S i z e [ Z ] + kk ] == v a l u e ) { 71 r e t u r n t r u e; 72 } 73 } 74 } 75 } 76 r e t u r n f a l s e ; 77 } 78 } ; 79 80 v o i d D i s t r i b u t o r : : e x e c u t e ( ) { 81 l o g g e r ( ) << ” R e c e i v i n g domain . . . ” << s t d : : e n d l ; 82 83 a u t o domainObs = i n s t a n c e . r e c e i v e (” domainIn ”) ; 84 t c u r = domainObs . timestamp ( ) ; 85 l o g g e r ( ) << ”Domain r e c e i v e d f o r t i m e ” << t c u r << s t d : : e n d l ; 86 87 s i z e t d o S i z e = domainObs . d a t a ( ) . s i z e ( ) ; 88 b o x S i z e = domainObs . d a t a ( ) . s h a p e ( ) ; 89 l o g g e r ( ) << ”Domain s i z e ” << d o S i z e << ” e l e m e n t s , ” <<

(38)

90 b o x S i z e [ 0 ] << ” by ” << b o x S i z e [ 1 ] << ” by ” <<

91 b o x S i z e [ 2 ] << s t d : : e n d l ;

92

93 a u t o domainObsPtr = domainObs . d a t a ( ) . e l e m e n t s ( ) ;

94 aggregateDomain . a s s i g n ( domainObsPtr , domainObsPtr + d o S i z e ) ;

95 l o g g e r ( ) << ” A g g r e g a t e domain a s s i g n e d ” << s t d : : e n d l ; 96 97 c o m p u t e A n d S e n d I n c r e m e n t a l B o u n d a r i e s ( aggregateDomain , 98 p r e v i o u s S o l i d ) ; 99 aggregateDomain . c l e a r ( ) ; 100 l o g g e r ( ) << ” D i s t r i b u t o r done ” << s t d : : e n d l ; 101 } 102 103 v o i d D i s t r i b u t o r : : c o m p u t e A n d S e n d I n c r e m e n t a l B o u n d a r i e s ( 104 c o n s t s t d : : v e c t o r & domain ,

105 s t d : : v e c t o r & p r e v i o u s S o l i d ) { 106 s i z e t l e n = domain . s i z e ( ) ; 107 // 0 − f l u i d , 1 − s o l i d , 2 − s u r f a c e 108 s t d : : v e c t o r i s S o l i d ( l e n , 0 ) ; 109 110 s i z e t n = 0 ; 111 f o r ( s i z e t i = 0 ; i < b o x S i z e [X ] ; i ++) { 112 f o r ( s i z e t j = 0 ; j < b o x S i z e [Y ] ; j ++) { 113 f o r ( s i z e t k = 0 ; k < b o x S i z e [ Z ] ; k++) { 114 i f( domain [ i ∗ b o x S i z e [Y] ∗ b o x S i z e [ Z ] + 115 j ∗ b o x S i z e [ Z ] + k ] != FLUID) { 116 i s S o l i d [ n ] = 2 ; 117 } 118 n++; 119 } 120 } 121 } 122 l o g g e r ( ) << ” i s S o l i d a s s i g n e d ” << s t d : : e n d l ; 123 124 b o o l e q u a l = t r u e; 125 i f ( i s S o l i d . s i z e ( ) == p r e v i o u s S o l i d . s i z e ( ) ) { 126 f o r( s i z e t l = 0 ; l < i s S o l i d . s i z e ( ) ; ++l ) { 127 i f ( i s S o l i d [ l ] != p r e v i o u s S o l i d [ l ] ) { 128 e q u a l = f a l s e; 129 b r e a k; 130 } 131 } 132 } e l s e { 133 e q u a l = f a l s e; 134 } 135 136 i f ( e q u a l ) {

137 // Send empty message t o s i g n a l : no change o r f u r t h e r work . 138 i s S o l i d . c l e a r ( ) ; 139 l o g g e r ( ) << ” S e n d i n g empty message : n o t h i n g t o do ” << 140 s t d : : e n d l ; 141 } e l s e { 142 // S t o r e t h e changed s o l i d mask 143 p r e v i o u s S o l i d = i s S o l i d ; 144 l o g g e r ( ) << ” S t o r i n g changed geometry ” << s t d : : e n d l ; 145 }

(39)

146 a u t o s o l i d r e s u l t = Data : : g r i d ( i s S o l i d . d a t a ( ) , { b o x S i z e [X] , 147 b o x S i z e [Y] , b o x S i z e [ Z ] } , {”x”, ”y”, ” z ”} ) ; 148 149 i n s t a n c e . s e n d (” geometryOut1 ”, Message ( t c u r , s o l i d r e s u l t ) ) ; 150 i n s t a n c e . s e n d (” geometryOut2 ”, Message ( t c u r , s o l i d r e s u l t ) ) ; 151 l o g g e r ( ) << ” D i s t r i b u t o r : d a t a s e n t ” << s t d : : e n d l ; 152 } 153 154 v o i d r u n D i s t r i b u t o r (i n t a r g c , c h a r ∗ a r g v [ ] ) { 155 I n s t a n c e i n s t a n c e ( a r g c , argv , { 156 { Op er at or : : F INIT , {” domainIn ”} } , 157 { Op er at or : : O F , {” geometryOut1 ”, 158 ” geometryOut2 ”}} 159 } ) ; 160 161 D i s t r i b u t o r d i s t r i b u t o r ( i n s t a n c e ) ; 162 163 w h i l e ( i n s t a n c e . r e u s e i n s t a n c e ( ) ) { 164 d i s t r i b u t o r . e x e c u t e ( ) ; 165 } 166 } 167 168 i n t main (i n t a r g c , c h a r ∗ a r g v [ ] ) { 169 r u n D i s t r i b u t o r ( a r g c , a r g v ) ; 170 r e t u r n EXIT SUCCESS ; 171 }

Adapting submodules of the in-stent restenosis model to enable uncertainty quantification using MUSCLE 3

Bachelor Informatica