Towards a Model-Based Experimental Design of the Maturation Process of Biohybrid Heart Valves⁎

Available online at www.sciencedirect.com

ScienceDirect IFAC PapersOnLine 52-26 (2019) 113–120

Towards a Model-Based Experimental Towards a Model-Based Experimental Towards a Model-Based Experimental Towards a Model-Based Experimental Design of the Maturation Process of Towards a Model-Based Experimental Design of the Maturation Process of Design of the Maturation Process of  Design of the Maturation Process of  Biohybrid Heart Valves Design of the Maturation Process of  Biohybrid Heart Valves  Biohybrid Heart Valves Biohybrid Heart Heart Valves Valves  Biohybrid Kirsten Voß ∗∗ Lorenz Pyta ∗∗ Jonas Gesenhues ∗∗ Petra Mela ∗∗ ∗∗

∗ Jonas Gesenhues ∗ Petra Mela ∗∗ Kirsten Voß ∗∗∗ Lorenz Pyta ∗ ∗ ∗ ∗∗ ∗∗ ∗ Petra Kirsten LorenzSchmitz-Rode Pyta Mela ∗∗ Dirk Abel ∗ Kirsten Voß VoßThomas Pyta ∗∗ Jonas Jonas Gesenhues Gesenhues Mela ∗∗ Thomas ∗ LorenzSchmitz-Rode ∗ Petra ∗∗ ∗∗ Dirk Abel ∗ Kirsten VoßThomas Lorenz Pyta Jonas Gesenhues Mela ∗∗ Dirk AbelPetra ∗ Schmitz-Rode Thomas Schmitz-Rode ∗∗ Dirk Abel ∗ ∗ Thomas Schmitz-Rode Dirk Abel ∗ Institute of Automatic Control, RWTH Aachen University, Germany ∗ Institute of Automatic Control, RWTH Aachen University, Germany ∗ Institute of Automatic RWTH (e-mail:Control, [email protected]) Control, RWTH Aachen Aachen University, University, Germany Germany (e-mail: [email protected]) ∗ Institute of Automatic ∗∗ of Automatic Control, RWTH Aachen University, Germany (e-mail: [email protected]) Institute of Applied Medical Engineering, RWTH Aachen University, ∗∗Institute (e-mail: [email protected]) Institute of Applied Medical Engineering, RWTH Aachen University, ∗∗ (e-mail: [email protected]) ∗∗ Institute of Applied Medical Engineering, RWTH Aachen University, Germany Engineering, RWTH Aachen University, ∗∗ Institute of Applied Medical Germany Institute of Applied Medical Germany Engineering, RWTH Aachen University, Germany Germany Abstract: Model-based experimental designs minimize the experimental effort while maximizAbstract: Model-based experimental designs minimize the experimental effort while maximizAbstract: Model-based experimental designs minimize experimental while ing the amount of analyzable experimental data. In thisthe paper, an initial effort mathematical model Abstract: Model-based experimental designs minimize the experimental effort while maximizmaximizing the amount of analyzable experimental data. In this paper, an initial mathematical model Abstract: Model-based experimental designs minimize the experimental effort while maximizing the amount of analyzable experimental data. In this paper, an initial mathematical model of the spatially distributed maturation process of tissue engineered heart valves is developed ing thespatially amount of analyzable experimental data.ofIntissue this paper, an initial model of the the distributed maturation process engineered heartmathematical valves is is developed ing amount of analyzable experimental data. this paper, anstate initial mathematical model of spatially distributed maturation of engineered heart valves forthe athemodel-based experimental design. process The model contains the variables amount, of spatially distributed maturation process ofIntissue tissue engineered heart valves cell is developed developed contains the state variables cell amount, for a model-based experimental design. The model of the spatially distributed maturation process of tissue engineered heart valves is developed for a model-based experimental design. The model contains the state variables cell amount, collagen concentration and elastin concentration. Based on the developed model,cell a variancefor a model-based experimental design. The model contains thedeveloped state variables amount, collagen concentration and elastin concentration. Based on the model, aa variancefor a model-based experimental design. The model contains thedeveloped state indicate variables cell amount, collagen concentration and elastin concentration. Based on the model, variancebased sensitivity analysis using Sobol’s method is performed. The results that all model collagen concentration and elastin concentration. Based on the developed model, a variancebased sensitivity analysis using Sobol’s method is performed. The results indicate that all model collagen concentration and elastin concentration. Based on the developed model, a variancebased sensitivity analysis using Sobol’s method is performed. The results indicate that all model parameters influence the model solution, while the variance of the parameters related to cell based sensitivity analysis method performed. results indicate related that all to model parameters influence the using modelSobol’s solution, whileis the varianceThe of the the parameters cell based sensitivity analysis using Sobol’s method isthe performed. The results indicate that all should model parameters influence the model solution, while the variance of parameters related to cell diffusion comprises anthe exceeding influence. Consequently, corresponding experiments parameters influence model solution, while variance of the parameters related to cell diffusion comprises an exceeding influence. Consequently, corresponding experiments should parameters influence model solution, while the variancecorresponding of the parameters related should to cell diffusion an influence. Consequently, experiments especiallycomprises focus on those parameters. diffusion comprises antheexceeding exceeding influence. Consequently, corresponding experiments should especially focus on those parameters. diffusion comprises an exceeding influence. Consequently, corresponding experiments should especially focus on those parameters. especially focus on those parameters. © 2019, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. especially on those parameters. Keywords: focus Distributed parameter systems, Tissue engineering, Numerical analysis, Keywords: Distributed parameter systems, Tissue engineering, Numerical analysis, Keywords: Distributed parameter systems, Tissue engineering, Numerical analysis, Mathematical models, Model-based experimental design, Sensitivity analysis Keywords: Distributed parameter systems, Tissue design, engineering, Numerical analysis, Mathematical models, Model-based experimental Sensitivity analysis Keywords: Distributed parameter systems, Tissue engineering, Numerical analysis, Mathematical models, Model-based experimental design, Sensitivity analysis Mathematical models, Model-based experimental design, Sensitivity analysis Mathematical models, Model-based experimental design, Sensitivity analysis to be spatially distributed, meaning that the tissue growth 1. INTRODUCTION 1. to distributed, meaning that 1. INTRODUCTION INTRODUCTION to be be spatially spatially distributed, meaning that the the tissue tissue growth growth properties vary across the heart valve. 1. INTRODUCTION to be spatially distributed, meaning that the tissue growth properties vary across the heart valve. INTRODUCTION be spatially distributed, meaning that the tissue growth properties vary across the heart valve. Tissue engineered1.implants are a valuable contribution to properties vary across the heart valve. For a robust implantable biohybrid heart Tissue engineered implants are valuable contribution properties varyproduction across the of heart valve. For production of biohybrid heart Tissue engineered implants implants are aaa valuable valuable contribution in the replacement of non-functional tissues and organs. Tissue engineered are contribution For aaa robust robust production of implantable implantable biohybrid heart valves it is necessary to standardize the complex matutissues and organs. in the replacement of non-functional For robust production of implantable biohybrid heart Tissue engineered implants are a valuable contribution valves it is necessary to standardize the complex matuin the replacement of non-functional tissues and organs. So-called biohybrid implants consist of a textile scaffold, a robust production of implantable biohybrid heart in the replacement of non-functional tissues andscaffold, organs. For valves it is necessary to standardize the complex maturation process while considering the uniqueness of every So-called biohybrid implants consist of a textile valves it is necessary to standardize the complex matuin the isreplacement of non-functional and organs. ration process while considering the uniqueness of every So-called biohybrid implants consistconnective oftissues textile scaffold, which covered with in vitro grown tissue. Due valves it is necessary to standardize the complex matuSo-called biohybrid implants consist of aa textile scaffold, ration process while considering the uniqueness of every implant. Therefore, the automation of the maturation which is covered with in vitro grown connective tissue. Due ration process while the considering the uniqueness of every So-called biohybrid consist of acombine textile scaffold, Therefore, automation of the maturation which is components, covered withimplants in vitro vitro grown grown connective tissue. Due to theiris biohybrid implants mechanration process while the considering the uniqueness of every which covered with in connective tissue. Due implant. implant. Therefore, automation of the maturation process of biohybrid heart valves is required. The main to their components, biohybrid implants combine mechanimplant. Therefore, the automation of the maturation which is covered with in vitro grown connective tissue. Due process of biohybrid heart valves is required. The main to their components, biohybrid implants combine mechanicaltheir stability and biological compatibility (Sodhani et al., implant. Therefore, the automation of the maturation to components, biohybrid implants combine mechanprocess of biohybrid heart valves is required. The main challenge of the automated production of these artificial ical stability and biological compatibility (Sodhani et al., process of biohybrid heart valves is required. The main to their components, biohybrid implants combine mechanchallenge of the automated production of these artificial ical stability and biological compatibility (Sodhani et al., 2018). Among other applications they offer great potential process of biohybrid heart valves is required. The main ical stability and biological compatibility (Sodhani et al., challenge of the automated production of these artificial heart valves is the high degree of individuality of each liv2018). Among other applications they offer great potential challenge of the automated production of these artificial ical stability and biological compatibility (Sodhani et al., heart valves is the high degree of individuality of each liv2018). Among other applications they offer offer greatbiohybrid potential in the area of artificial heart valves, because challenge of the automated production of these artificial 2018). Among other applications they great potential heart valves is the high degree of individuality of each living tissue and therefore each growing heart valve through in the area of artificial heart valves, because biohybrid heart valves is therefore the high degree of individuality of each liv2018). Among other applications they offer great potential ing tissue and each growing heart valve through in the area of artificial heart valves, because biohybrid heart valves are able to withstand mechanical stress while heart valves is the high degree of individuality of each livin thevalves area are of artificial heart valves, becausestress biohybrid ing tissue and therefore each growing heart valve through e.g. individual speeds of growth. A model-based control heart able to withstand mechanical while ing tissue and therefore each growing heart valve through in the area of artificial heart valves, because biohybrid e.g. individual speeds of growth. A model-based control heart valves are able to withstand mechanical stress while being valves highly are biologically compatible. Furthermore, the liv- ing tissue and therefore each growing heart valve through heart able to withstand mechanical stress while e.g. individual speeds of growth. A model-based control scheme is aa possible to this problem: Depending being highly biologically compatible. Furthermore, the livindividual speedssolution of growth. A model-based control heart valves are able to the withstand stress while scheme is solution to this problem: Depending being highly biologically compatible. Furthermore, the living tissue which covers scaffoldmechanical isFurthermore, developed from the e.g. e.g. individual speeds of growth. A model-based control being highly biologically compatible. the livscheme is aa possible possible solution to this problem: Depending on whether the parameters vary during the maturation ing tissue which covers the scaffold is developed from the scheme is possible solution to this problem: Depending being highly biologically compatible. the livwhether the parameters vary during the maturation ing tissue which covers the scaffold isFurthermore, developed from the on patient’s own cells (Flanagan et al., 2007). For this reason, scheme is a possible solution to this problem: Depending ing tissue which covers the scaffold is developed from the on whether the parameters vary during the maturation process, the patient specific model parameters can be patient’s own cells (Flanagan et al., 2007). For this reason, whether the parameters vary during the maturation ing tissueown which covers theable scaffold is developed the on process, the patient specific model parameters can be patient’s own cells (Flanagan ettoal., al., 2007). For thisfrom reason, biohybrid heart valves are grow andFor regenerate and on whether the parameters vary during the maturation patient’s cells (Flanagan et 2007). this reason, process, the patient specific model parameters can be identified either offline ahead or online during the process. biohybrid heart valves are able to grow and regenerate and process, the patient specific model parameters can be patient’s own cells (Flanagan et al., 2007). For this reason, identified either offline ahead or online during the process. biohybrid heart valves are able to grow and regenerate and therefore adapt to their surrounding. Correspondingly, the patient specific model parameters can be biohybrid heart valves aresurrounding. able to grow Correspondingly, and regenerate andaa process, identified either offline ahead or online during the process. In doing so, the model-based control is able to adapt to therefore adapt to their identified either offline ahead or online during the process. biohybrid heart valves are able to grow and regenerate and In doing so, the model-based control is able to adapt to therefore adapt to their surrounding. Correspondingly, a biohybrid heart valve is a complex system, which reacts identified either offline ahead orcontrol online during the process. therefore adapt to their surrounding. Correspondingly, a the In doing so, the model-based is able to adapt to individual heart valve. Thus, mathematical models are biohybrid heart valve is a complex system, which reacts In doing so, the model-based control is able to adapt to therefore adapt to their surrounding. Correspondingly, a the individual heart valve. Thus, mathematical models are biohybrid heart valve is a complex system, which reacts to and interacts with its environment. Though there exist In doing so, the model-based control is able to adapt to biohybrid heart with valveits is environment. a complex system, which reacts the individual heart valve. Thus, mathematical models are a great contribution to the understanding of this maturato and interacts Though there exist the individual heart valve. Thus, mathematical models are biohybrid heart valve is a complex system, which reacts a great contribution to the understanding of this maturato andapproaches interacts with with itsasenvironment. environment. Though there exist other suchits in situ tissue engineering (e.g. the individual heart valve. mathematical models are to and interacts Though there exist a great great contribution to theThus, understanding of this this maturation process and the development of corresponding robust other approaches such as in situ tissue engineering (e.g. a contribution to the understanding of maturato andapproaches interacts with Though there process and the to development of corresponding robust other approaches suchitsthe asenvironment. inmaturation situ tissue tissue engineering (e.g. Bouten et al., 2018), process of a exist bio- tion ation great contribution the understanding of thisamaturaother such as in situ engineering process and development of robust control algorithms. the validation of such model, Bouten et al., 2018), the maturation process of aa (e.g. biotion process and the the For development of corresponding corresponding robust other approaches such as in situ tissue engineering (e.g. control algorithms. For the validation of such a model, Bouten et al., 2018), the maturation process of biohybrid heart valve is usually executed in a bioreactor. tion process and the development of corresponding robust Bouten et al., 2018), the maturation process of a biocontrol algorithms. For the validation of such a model, experimental data are required, but experiments are time hybrid heart valve is usually executed in a bioreactor. control algorithms. For the validation of such a model, Bouten et al.,valve 2018), process of aetbiodata are required, but experiments are time hybrid heart valve is the usually executed in bioreactor. One potential bioreactor is maturation described by Flanagan al. experimental control algorithms. For the validation of such a model, hybrid heart is usually executed in aa bioreactor. experimental data are required, but experiments are consuming and expensive: The maturation process of a One potential bioreactor is described by Flanagan et al. data are required, but experiments are time time hybrid heart is usually executed consuming and expensive: The maturation process of One potential bioreactor is described described by in Flanagan et al. al. experimental (2007): While valve the heart valve is perfused bya abioreactor. pulsating experimental data are required, but experiments are time One potential bioreactor is by Flanagan et consumingheart and valve expensive: The maturation process of a a biohybrid has a duration of three weeks (Mor(2007): While the heart valve is perfused by a pulsating consuming and expensive: The maturation process of a One bioreactor is to described by et al. biohybrid heart has aa The duration of three weeks (Mor(2007): While the heart valve valve is perfused perfused by pulsating flow potential of While nutrition solution imitate theFlanagan cardiac cycle, consuming and valve expensive: maturation process of a (2007): the heart is by aa pulsating biohybrid heart valve has duration of three weeks (Moreira et al., 2016). Apart from this, due to the vulnerability flow of nutrition solution to imitate the cardiac cycle, biohybrid heart valve has a duration of three weeks (Mor(2007): thebe heart valve is perfused a actuating pulsating et al., heart 2016). Apart from this, dueofto the vulnerability flow of While nutrition solution to imitate imitate the by cardiac cycle, eira the process can manipulated by multiple biohybrid has acontamination, duration weeks (Morflow of nutrition solution to the cardiac cycle, eira et 2016). Apart from this, the vulnerability of the growing tissue to the maturation actuating the process can be manipulated by multiple eira et al., al., 2016).valve Apart from this, due due to tothree the vulnerability flow of nutrition solution to imitate the cardiac cycle, of the growing tissue to contamination, the maturation the process can be manipulated by multiple actuating variables. Mechanical stimuli such as hydrostatic pressure eira et al., 2016). Apart from this, due to the vulnerability the process can be manipulated byhydrostatic multiple actuating of the growing tissue contamination, the maturation process frequently fails to finish successfully. Consequently, variables. Mechanical stimuli such as pressure the growing tissue to finish contamination, the maturation the can bystimuli multiple actuating process frequently fails to successfully. Consequently, variables. Mechanical stimuli such as as hydrostatic pressure and process flow rate as be wellmanipulated as biological through the of of thea growing tissue to contamination, the maturation variables. Mechanical stimuli such hydrostatic pressure process frequently fails finish successfully. Consequently, only low number of accomplished experiments can be and flow rate as well as biological stimuli through the process frequently fails to finish successfully. Consequently, variables. Mechanical stimuli such as hydrostatic pressure only a low number of accomplished experiments can be and flow rate as well as biological stimuli through the composition of the nutrient solution can be applied. As an process frequently fails to finish successfully. Consequently, and flow rate as well as biological stimuli through the only a low number of accomplished experiments can be performed. In addition, the only present possibility of be applied. As an composition of the nutrient solution can only a low number of accomplished experiments can bea and flow rate as well as biological stimuli through the performed. In addition, the only present possibility of composition of the nutrient solution can be applied. As an additional complexity, the maturation process is expected a low number of accomplished experiments can bea composition of the nutrient solution canprocess be applied. As an only performed. In addition, the only present possibility of a state analysis is a histological analysis of the heart valve, additional complexity, the maturation is expected performed. In addition, the only present possibility of a composition of the nutrient solution can be applied. As an state analysis is a histological analysis of the heart valve, additional complexity, the maturation process is expected performed. In addition, the only present possibility of a additional complexity, the maturation process is expected which state analysis is a histological analysis of the heart valve, destroys it. Therefore, such a snap-shot of the matustate analysis is a histological analysis of the heart valve,  This workcomplexity, additional the maturation process is expected which destroys it. Therefore, such a snap-shot of the matufunded by the Deutsche Forschungsgemeinschaft state analysis is a histological analysis of the heart valve,  This work was which destroys it. Therefore, such a snap-shot of the maturating heart valve results in aasuch low of experimental was funded by the Deutsche Forschungsgemeinschaft which destroys it. Therefore, aamount snap-shot of the matu This German rating heart valve results in low of experimental (DFG, Research – 403043858. work funded by Forschungsgemeinschaft  which destroys it. Therefore, aamount snap-shot of the matuThis German work was was funded Foundation) by the the Deutsche Deutsche Forschungsgemeinschaft rating heart valve results in aasuch low amount of experimental (DFG, Research Foundation) – 403043858.  rating heart valve results in low amount of experimental (DFG, Research – This German work was funded Foundation) by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 403043858. 403043858. rating heart valve results in a low amount of experimental

(DFG, German Research Foundation) –Federation 403043858. 2405-8963 © 2019, IFAC (International of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Peer review under responsibility of International Federation of Automatic Control. 10.1016/j.ifacol.2019.12.245

114

Kirsten Voß et al. / IFAC PapersOnLine 52-26 (2019) 113–120

data. For the mentioned reasons, only a low number of experiments which result in a small amount of analyzable data can be performed. To obtain the maximum benefit from a minimum number of experiments, it is important to perform a wellstructured design of experiments (Franceschini and Macchietto, 2008). To gain further knowledge of the process and the most influential model parameters, a sensitivity analysis of an initial model of the process can be performed and provide a basis for a systematic experimental design. In this paper, a simplified mathematical model of the in vitro maturation process of biohybrid heart valves is developed. The one-dimensional model consists of a coupled system of partial and ordinary differential equations to represent cell proliferation dynamics as well as the development of the concentrations of the proteins collagen and elastin to capture the most important dynamics during the maturation process. To obtain the influence of the model parameters regarding the state variables cell amount and protein concentrations, a corresponding sensitivity analysis using Sobol’s method is performed. This paper is structured as follows: In Section 2 the basic concepts of the modeled phenomenons are introduced. Furthermore, relevant existing modeling approaches are presented. In Section 3 a potential spatially distributed mathematical model is developed. The implementation and a numerical solution of the model is described in Section 4. In Section 5 the applied sensitivity analysis method is described. The results of the sensitivity analysis are presented in Section 6, while the discussion of the results and conclusions for the design of experiments follow in Section 7. A summary concludes the paper in Section 8. 2. STATE OF THE ART For the development of a process model, it is essential to understand the process itself. During the maturation process of biohybrid implants, connective tissue, which consists of valvular interstitial cells (VICs) and an extracellular matrix (ECM), is formed: The initial cells which are seeded on the scaffold proliferate and simultaneously build the ECM through biosynthesis of ECM components. VICs hold characteristics of the phenotypes of fibroblasts and smooth muscle cells, whereas it is not clear whether VICs belong to the phenotype of myofibroblasts or whether fibroblasts and smooth muscle cells coexist in heart valve tissue (Mulholland and Gotlieb, 1996). Due to limitations of nutrition transports to the inner cells and waste removal out of the tissue, the thickness of the tissue in vivo and in vitro is limited to 100 - 200 µm (Carmeliet and Jain, 2000). The ECM consists of interstitial fluid, proteins and polysaccharides (Alberts et al., 2015). In human heart valve tissue, the proteins collagen and elastin play an essential role: Collagen and elastin are, respectively, responsible for the stress-bearing and the elasticity of a heart valve. Through a three-layered architecture with increased elastin concentration on the inflow side and increased collagen concentration on the outflow side, a heart valve is able to withstand the dynamics of the cardiac cycle (Schoen, 2012). Because collagen and elastin are expected to be equally important for the robustness of biohybrid

heart valves, both of these proteins have to be considered in a mathematical model besides cells. To the best of the authors’ knowledge, models of the maturation process of biohybrid implants which consider cells as well as the mentioned ECM components do not exist yet. However, modeling approaches with similar purposes can be found in the literature: Hossain et al. (2015) model the chondrocyte cell growth and cell-nutrition interaction in a perfusion bioreactor. A model of cell population growth in a one dimensional scaffold is analyzed by Shakeel (2013). Bouten et al. (2018) develop a mechanobiological model of cell and collagen dynamics for the development of in situ tissue engineered heart valves. Wound healing models (e.g. Olsen et al., 1995) are developed to predict the tissue reconstruction of wounded tissue. The dominant general approach for modeling cell population growth is the one-dimensional Fisher–Kolmogorov equation (FKE), which is a non-linear diffusion-reactionequation. Through space and transport limitations, the cell proliferation can be modeled as a saturation growth curve. Besides the proliferation, cell diffusion is an important phenomenon during tissue formation, because the cells need to move to create space for new cells. The FKE can be written as   ∂n ∂ ∂n = D + rn (1 − n) , (1) ∂t ∂x ∂x

where n is the population size, D is the diffusion coefficient of the observed population and r is the corresponding growth rate (Fisher, 1937; Kolmogoroff et al., 1988). The independent variables t and x indicate that the population size is distributed in time t and in one spatial dimension x. The equation is nondimensionalized such that a value of n = 1 corresponds to the steady state population size, e.g. the normal cell amount in connective tissue.

For the development of a model of the maturation process of biohybrid heart valves, models of similar processes are combined to address the present application. For the present application especially wound healing models are relevant, because the appearing phenomenons of wound healing processes are similar to the ones of the growth of biohybrid implants. During the wound healing process, new tissue is formed while cell and protein dynamics are influenced by stimuli of the surrounding healthy tissue. Additionally, in vitro VIC proliferation holds the same characteristics as fibroblasts proliferation during wound healing (Benton et al., 2009). A popular model of tissue growth in wound healing which is based on the FKE was developed by Olsen et al. (1995, 1997). 3. MATHEMATICAL MODEL Taking into account the lack of knowledge of the specific subphenotypes of VICs in the literature, a general VIC amount nV is chosen as a state variable for the process model. Based on the previous considerations on the threelayered architecture of a heart valve, collagen concentration cC and elastin concentration cE are chosen as additional state variables. In Fig. 1 a simplified scaffold, which consists of three leaflets like a human heart valve does, and the coordinate system of the model are presented. Even though in reality the tissue grows in every direction



Kirsten Voß et al. / IFAC PapersOnLine 52-26 (2019) 113–120

orthogonal to the scaffold, for a first modeling approach only one direction is considered. The spatial axis x of the model represents the distance between the scaffold and the regarded position, while the axis of the state variables represents a non-spatial dimension. The spatial axis x is sufficient for an initial model: As cell cultures aim to minimize their surface tension (Albert and Schwarz, 2014, 2016), the tissue growth along the scaffold can be assumed to be uniform if only the inner area of a leaflet is considered. Consequently, there is no spatial gradient of the VIC amount in the direction of the y and z axis. However, the model can easily be expanded to additional dimensions if required.

y y z z x

nV , cC , cE

x

∂ ∂nV = ∂t ∂x



115

DV · e

+ rV nV



γ



1−

cE cE,0

xmax



nV nV,max

 ∂nV ∂x  − d V nV .

(4)

Here, rV is the growth rate, dV is the decay constant, nV,max is the maximum VIC amount and cE,0 is the steady state elastin concentration. The diffusion coefficient of the VICs DV is weakened by the exponential term to model the restricted diffusion. The diffusion coefficient γ defines the strength of this effect. A functional diagram of the model is presented in Fig. 2. As the proteins are synthesized by the VICs, the VIC amount influences the collagen and elastin concentration. In turn, elastin as a representative of the ECM influences the diffusion of the VICs and therefore the VIC amount, while the collagen concentration is a pure output variable. nV

Scaffold

−1

cC

x cE

Fig. 1. Coordinate system of the one-dimensional model of the maturation process of biohybrid heart valves.

Fig. 2. Interaction of the state variables of the model: the VIC amount nV , the collagen concentration cC and the elastin concentration cE .

The following model is based on a model of wound healing, which is developed by Olsen et al. (1995, 1997). Olsen’s model is adjusted for the present application as follows: During the maturation process, there is no displacement of the ECM and therefore also no convection. Furthermore, besides cells and collagen, elastin is considered as state variable, because this protein is important for the architecture of heart valves as well. The dynamics of elastin are the same as the dynamics of collagen. The model equations of the proteins result as follows:   2 ∂cC rC RC = (2) 2 + c2 − dC cC nV , ∂t RC C

To find a suitable solution from the solution space of the system of equations (2) to (4), initial conditions for all equations are required. In addition, two boundary conditions have to be defined for the solution of the second order partial equation of the VICs (4). The boundary and initial conditions are determined as follows: Due to the symmetric growth of the tissue in every direction orthogonal to the scaffold, a Neumann boundary condition for the VIC equation (4) at x = 0 is chosen:  ∂nV  = 0. (5) ∂x 

∂cE = ∂t





2 rE RE 2 + c2 − dE cE nV . RE E

(3)

The development of collagen and elastin is determined by the respective growth rates rC and rE , the saturation constants RC and RE and the decay constants dC and dE . To the best of the authors’ knowledge, little is known about the parameters of VIC proliferation such as the growth rate. Because the corresponding cellular phenotype is fibroblast-like, dynamics and parameter values of fibroblast proliferation from Olsen’s wound healing model are used to model the development of the VIC amount. In contrast to cells present in wounds, the VICs of the maturation process can not diffuse freely. Instead, diffusion is greater at the point where ECM (here represented by the elastin concentration) is already built. Similar to the nonlinear diffusion model of Shakeel (2013), an exponential term is used to model this phenomenon, because this term fulfills the mentioned requirements. The development of the VIC amount results as follows:

x=0

Furthermore, due to the limited thickness of a grown tissue, represented by the spatial boundary xmax , the following Dirichlet boundary condition is chosen: (6) nV (x = xmax ) = 0. At the beginning of the maturation process, the scaffold is covered with cells. Considering that the size of the seeded cells is similar to the size of fibroblasts, which hold an approximate diameter of 16 µm (Angello et al., 1987), the initial thickness of the cell layer covering the scaffold is set to 22.5 µm. The initial condition of the VICs can therefore be expressed through a modified Heaviside step function H(·) multiplied by the steady state VIC amount nV,0 , which equals the initial amount: (7) nV (x, t = 0) = nV,0 (1 − H (x − 22.5 µm)) . ECM and therefore proteins are not present at the ginning of the maturation process (see Flanagan et 2007). The initial values of the protein concentrations therefore set to zero: cC (x, t = 0) = 0, cE (x, t = 0) = 0.

beal., are (8) (9)

Kirsten Voß et al. / IFAC PapersOnLine 52-26 (2019) 113–120

Table 1. Parameters of the presented model and values used in the simulation derived from the literature.

4. IMPLEMENTATION AND SIMULATION The corresponding parameter values used for the simulation are determined by Olsen et al. (1995), Bashey et al. (1967), Shakeel (2013), Sherratt (2009) and Carmeliet and Jain (2000). The parameter values of the elastin equation (3) are calculated in an analogous manner to the calculation of the parameters of the collagen equation (2) performed by Olsen et al. (1995). Based on preliminary considerations on the delimitation of the growing tissue and the nutrient solution, the diffusion parameter γ in this application is expected to be higher compared to Shakeel (2013). Based on a qualitative sensitivity analysis of the diffusion parameter, γ = 5 is chosen for the simulation. The parameter values used for the simulation are presented in table 1. The model is implemented in MATLAB (2018). The spatial derivative comprised in the VIC equation (4) is implemented using the method of Crank and Nicolson (1996) for spatial discretization, while the time discretization is proceeded automatically by the solver. For time discretization the 4th order Runge-Kutta method is chosen, because this explicit method comprises less computational effort than implicit methods and is more accurate compared to other explicit methods such as Euler’s method (Abel and Bollig, 2006). To consider all time constants of the model, a sampling time of 3 minutes is chosen. Furthermore a distance of 1 µm is used for the spatial discretization points to obtain a sufficient resolution. For simplification and to transform all computed quantities to the same order, the model is nondimensionalized such that a value of 1 of the state variables, respectively, equals the normal cell amount or protein concentration in human connective tissue. Furthermore, the model is scaled to the maximum possible tissue thickness of xmax = 150 µm and the simulated maturation process duration of three weeks. The non-differentiable initial condition of the

0.8 0.6 0.4

t

0.2 0

0

0.2 0.4 0.6 0.8 1 Distance to scaffold x

00

0.5 1 Distance to scaffold x

21 14 7 Time t 0 [ d]

Collagen concentration cC

VIC amount nV

0.5

The development of the spatial distributions of the VIC amount and the collagen and elastin concentrations during the maturation process resulting from the simulation are presented in Fig. 3. Due to the chosen boundary conditions, the solutions are similar to wave front solutions. However, in contrast to real wave fronts, due to the nonconstant diffusion term, the solutions do not show a constant front velocity. According to the model solution, the VICs immediately diffuse into the free space while synthesizing collagen and protein. In doing so, the cells slightly diffuse into free space without already existing ECM (compare e.g. the distributions after three days). Conceivably, cells move to the edge of the tissue to synthesize ECM components and therefore tissue growth. Because experimental data are not available at this time, the simulation results can not be validated. However, the simulation results do not show any indication of significant

0.8

0

t 0

0.2 0.4 0.6 0.8 1 Distance to scaffold x ), 3 days (

), 6 days (

·10−2 1.6 1.2 0.8 0.4 00

5 1 · 10−4 g mm−3 3 · 10−5 g mm−3 1.7 · 10−8 mm2 s−1 7.5961 · 10−5 d−1 mm3 6.3301 · 10−6 d−1 mm3 0.0231 d−1 10 mm−3 1000 mm−3 3 · 10−4 g mm−3 9 · 10−5 g mm−3 8.4401 · 10−9 g d−1 2.1100 · 10−10 g d−1 0.0233 d−1 150 µm

VIC distribution (7) is approximated by a sigmoid function to obtain dual differentiability.

1.2

0.4

Value

Diffusion parameter Steady state collagen concentration Steady state elastin concentration Diffusion coefficient Decay constant of collagen per cell Decay constant of elastin per cell Decay constant of VICs Steady state VIC amount Maximum VIC amount Saturation constant of collagen Saturation constant of elastin Growth rate of collagen Growth rate of elastin Growth rate of VICs Maximum tissue thickness

−2 1.6 ·10

(a) Spatial distributions of the state variables after 0 days ( ( ) and 21 days ( ).

1

γ cC,0 cE,0 DV dC dE dV nV,0 nV,max RC RE rC rE rV xmax

Elastin concentration cE

Collagen concentration cC

VIC amount nV

1

Symbol Description

0.5 1 Distance to scaffold x

21 14 7 Time t 0 [ d]

), 9 days (

Elastin concentration cE

116

−3 1.5 ·10

1 0.5 0

t 0

0.2 0.4 0.6 0.8 1 Distance to scaffold x

), 12 days (

), 15 days (

), 18 days

·10−3 1.5 1 0.5 00

0.5 1 Distance to scaffold x

(b) Surface plots of the progression of the state variables during three weeks.

Fig. 3. Progression of the spatially distributed state variables during a period of three weeks.

21 14 7 Time t 0 [ d]

Kirsten Voß et al. / IFAC PapersOnLine 52-26 (2019) 113–120

modeling errors. Therefore, the resulting progression of the state variables could represent the tissue growth in a realistic manner. Overall, the simulation results demonstrate that the developed model is a suitable first attempt for an initial analysis and a model-based experimental design. The computing time for the simulation on a standard laptop is approximately 2 s. Considering the large time constants of the model, the developed model simulation is already suitable for a model-based control of the maturation process regarding its computational effort. 5. SENSITIVITY ANALYSIS The aim of the sensitivity analysis is to determine the connection between the model parameters and the solution of the mathematical model. Results of the analysis allow conclusions to be drawn about parameters on which the experiments need to be focused on and parameters which can be treated as constant due to their insignificant influence on the model solution. For this reason, the responsiveness of the model to parameter variations and uncertainties have to be analyzed (Cukier et al., 1978). For the global sensitivity analysis of non-linear systems, wellestablished approaches are the Fourier-Amplitude Sensitivity Test (FAST) and the method of Sobol’ (Henkel et al., 2012). The latter requires more computational effort, but explicit sensitivity indices can be calculated. Furthermore, this model-independent sensitivity method is known to be robust and accurate (Henkel et al., 2012). Due to this reasons, Sobol’s method is used to analyze the influence of parameter variation and uncertainties on the solution of the presented model. The sensitivity analysis is performed for each of the three state variables at each discrete spatial position and at each point in time. The grid points at t = 0 and x = xmax are not considered, because the values at these positions are specified by the initial conditions and the Dirichlet boundary condition. Consequently, for each state variable with J spatial positions and K discrete points of time a matrix (J − 1) × (K − 1) of sensitivity indices is received. During the sampling for the sensitivity analysis, all model parameters are varied except from the growth rates rV , rC and rE as well as the maximum VIC amount nV,max . The growth rates are each determined by the steady state values and the decay rates (see Olsen et al., 1995), while the maximum VIC amount is specified by spatial limitation. To attenuate the variation of the exponential term, instead of γ the whole exponential term γˆ = eγ is varied. As the real variation of the model parameters is not known to this point, the used parameter spaces equate to the values taken from the literature ±50 %. To determine whether a parameter value can be treated as constant, e.g. to a literature value, the total Sobol index ST i is the crucial one (Zhang et al., 2015). Therefore, this sensitivity index is determined and analyzed. The sensitivity analysis is implemented in MATLAB (2018) using the open-source SAFE-Toolbox (Pianosi et al., 2015). The functions of the toolbox are adapted to the parallel calculation of the sensitivity indices of twodimensional output variables. To determine the required number of model evaluations, a convergence analysis of the sensitivity indices is performed. Finally, the development of the sensitivity of each state variable on the model parameters over time is analyzed.

117

6. RESULTS Exemplary results of the convergence analysis of the sensitivity indices are presented in Fig. 4. Illustrated is the development of the total sensitivity indices at one spatial and temporal grid point regarding an increasing sample size. The convergence analyses of the other grid points yield to the same results. All total sensitivity indices converge at an approximate sample size of 6000. Therefore, this sample size is used for the determination of the sensitivity indices. Consequently, for the analysis of the sensitivity of the state variables regarding the corresponding model parameters, 60 000 model evaluations are performed (see Saltelli et al. (2010) for details). 0.8 0.6 ST i,V (Xi )



0.4 0.2 0

−0.2

1000

2000

3000 4000 Sample size

5000

6000

Fig. 4. Exemplary results of the convergence analysis of the total sensitivity indices of the analyzed parameters γˆ ( ), cC,0 ( ), cE,0 ( ), DV ( ), dC ( ), dE ( ), dV ( ), nV,0 ( ), RC ( ) and RE ( ) regarding the VIC amount at the grid point t = 0.5 and x = 0.5. The temporal development of the total sensitivity indices of the VIC amount ST i,V , of the collagen concentration ST i,C and of the elastin concentration ST i,E resulting from the sensitivity analysis is illustrated in Fig. 5. The spatial means of these sensitivity indices are presented, while spatial variations, meaning the lower and upper value along the spatial axis, are represented by colored areas. In addition, the values of the means as well as the spatial minima and maxima at the beginning and at the end of the simulated maturation process are presented in table 2. The highest mean total sensitivity indices (mSTis) regarding the VIC amount correspond to the adjusted diffusion parameter γˆ and the diffusion coefficient DV , see Fig. 5a. The mean total sensitivity index (mSTi) of γ decreases during the maturation process from 0.7628 to 0.4195. Similar to the mSTi of the adjusted diffusion parameter, the mSTi of the diffusion coefficient DV regarding the VIC amount decreases as well, in this case from 0.7537 to 0.4091. In contrast, the average total effect of the normal VIC amount nV,0 increases from 0.0897 to 0.1845. The mSTis of all other analyzed parameters regarding the VIC amount remain below 0.1. Just like the sensitivity of the VIC amount, the development of the collagen concentration mainly depends on the adjusted diffusion parameter γˆ and the diffusion coefficient Dn , see Fig. 5b. The mSTi of γˆ decreases from 0.8789 to 0.7045 while the mSTi of DV decreases from 0.8726 to 0.6377. The mSTi of the normal collagen concentration cC,0 slightly increases in the beginning of the maturation

Kirsten Voß et al. / IFAC PapersOnLine 52-26 (2019) 113–120

118

Having the highest impact on the elastin concentration as well, the mSTi of γˆ decreases from 0.9291 down to 0.7795, see Fig. 5c. The mSTi of DV decreases from 0.9014 down to 0.7689. Except from the mSTi of the elastin decay rate dE , all other mSTis increase from approximately 0.3 to 0.55. The mSTi of dE holds the lowest values while increasing from 0.2350 to 0.5405.

1

ST i,V (Xi )

0.8 0.6 0.4

Table 2. Total sensitivity indices of the analyzed parameters regarding the spatially distributed state variables at the beginning and at the end of the maturation process.

0.2 0

0

3

6

9 12 Time t [ d]

15

18

21

(a) Total sensitivity indices of the analyzed parameters regarding the VIC amount.

(a) Total sensitivity indices of the analyzed parameters regarding the VIC amount. 1

ST i,C (Xi )

0.8 0.6 0.4 0.2 0

0

3

6

9 12 Time t [ d]

15

18

21

(b) Total sensitivity indices of the analyzed parameters regarding the collagen concentration.

STi,V (Xi ) t=0d Xi

mean

min

max

mean

min

max

γ ˆ cC,0 cE,0 DV dC dE dV nV,0 RC RE

0.7628 0.0196 0.0196 0.7537 0.0196 0.0196 0.0197 0.0897 0.0196 0.0196

0 0 0 0 0 0 0 0 0 0

0.9999 0.0610 0.0610 0.9974 0.0610 0.0610 0.0612 0.5211 0.0610 0.0610

0.4195 0.0017 0.0017 0.4091 0.0017 0.0016 0.0013 0.1845 0.0017 0.0017

0.0042 0 0 0.0049 0 0 0 0.0435 0 0

0.7303 0.0101 0.0102 0.7132 0.0101 0.0101 0.0080 0.5171 0.0101 0.0101

(b) Total sensitivity indices of the analyzed parameters regarding the collagen equation. STi,C (Xi ) t=0d

1

ST i,E (Xi )

0.8 0.6 0.4 0.2 0

t = 21 d

0

3

6

9 12 Time t [ d]

15

18

21

(c) Total sensitivity indices of the analyzed parameters regarding the collagen concentration.

Fig. 5. Total sensitivity indices of the analyzed parameters γˆ ( ), cC,0 ( ), cE,0 ( ), DV ( ), dC ( ), dE ( ), dV ( ), nV,0 ( ), RC ( ), RE ( ), and the corresponding spatial lower and upper bounds. process from 0.4825 to a maximum of 0.5305 and decreases afterwards down to 0.4974. The dependency of the collagen concentration on the collagen decay rate dC slightly decreases with a corresponding mSTi of 0.4835 and an mSTi of 0.4586 in the end of the process. In contrast, the mSTi of the normal VIC amount nV,0 slightly increases from 0.4332 to 0.4741. All other mSTis regarding the collagen concentration decrease from approximately 0.46 down to 0.37.

t = 21 d

Xi

mean

min

max

mean

min

max

γ ˆ cC,0 cE,0 DV dC dE dV nV,0 RC RE

0.8789 0.4825 0.4551 0.8726 0.4835 0.4551 0.4551 0.4332 0.4538 0.4551

0.5109 0.1995 0.2915 0.5111 0.3439 0.2915 0.2915 0.2769 0.2761 0.2915

1.0000 0.7234 0.5323 0.9993 0.6499 0.5323 0.5324 0.6501 0.5467 0.5323

0.7045 0.4974 0.3704 0.6377 0.4586 0.3704 0.3704 0.4741 0.3761 0.3704

0.5104 0.3667 0.2904 0.5100 0.3529 0.2904 0.2901 0.3674 0.2733 0.2904

0.8997 0.7218 0.5099 0.7968 0.6471 0.5099 0.5101 0.6527 0.5456 0.5099

(c) Total sensitivity indices of the analyzed parameters regarding the elastin equation. STi,E (Xi ) t=0d

t = 21 d

Xi

mean

min

max

mean

min

max

γ ˆ cC,0 cE,0 DV dC dE dV nV,0 RC RE

0.9291 0.3192 0.2605 0.9014 0.3192 0.2350 0.3193 0.2957 0.3192 0.3141

0.7158 0.0682 0 0.7196 0.0682 0 0.0686 0.0236 0.0682 0.0295

0.9999 0.8459 0.8721 0.9976 0.8459 0.8529 0.8459 0.8558 0.8459 0.8502

0.7795 0.5571 0.5790 0.7689 0.5571 0.5405 0.5556 0.5631 0.5571 0.5668

0.7023 0.3825 0.4015 0.7130 0.3825 0.3292 0.3773 0.3839 0.3825 0.4036

0.8440 0.8442 0.8701 0.8442 0.8442 0.8511 0.8440 0.8543 0.8442 0.8485



Kirsten Voß et al. / IFAC PapersOnLine 52-26 (2019) 113–120

7. DISCUSSION First of all, considering the results of the sensitivity analysis, the influence of the parameters of the collagen equation (2) on the other two state variables, especially the elastin concentration, is questionable, because the collagen concentration is a pure output value, see Fig. 2. A reason for this phenomenon may be the neglected influence of the temporal history of the amount and concentrations of the state variables on the values at a certain point of time. In addition, the influence of the VIC amount of neighboring discrete locations is not considered in the sensitivity analysis. Those two simplifications could be the reason for an overestimation of the influence of the analyzed parameters. However, the overall ratio of the sensitivity indices can still be used to determine important parameters: All of the three state variables are especially sensitive to a variance or uncertainty of the diffusion coefficient Dn and the adjusted diffusion parameter γˆ and therefore γ. As a result, the corresponding initial experiments should put special emphasis on these two parameters, because the variances of these parameters highly influence the development of the VIC amount as well as the protein concentrations during the whole maturation process. For the exact determination of the two diffusion parameters, a microscopic analysis as performed by Bard and Hay (1975) can support those experiments. Besides Dn and γˆ , all other analyzed parameters influence the development of the state variables as well, meaning that none of the parameters can be neglected and simply set to the literature value. For example, the VIC amount and the collagen concentration each depend on, respectively, their normal amount and normal concentration. Experiments regarding the relationship of those steady state values in a patient’s connective tissue and the values in the biohybrid heart valve are beneficial: If a direct relationship exists, these parameter values can be determined offline ahead of the maturation process by analyzing a probe of the patient’s connective tissue. The transmission of wound healing models to the maturation process of biohybrid heart valves, especially the similarity of the dynamics of fibroblasts and VICs, needs to be evaluated. In addition, the potential necessity to consider different subphenotypes like fibroblasts and myofibroblasts in the process model instead of just one cell type has to be studied, because e.g. it is unclear whether all cell types, especially myofibroblasts, show random motion (Olsen et al., 1995) as represented by the diffusion term. To analyze the model parameters, it is advantageous to execute the maturation process of the biohybrid implants in a static nutrition solution meaning that there is no flow and therefore no mechanical influence such as convection. As the biohybrid heart valves do not open and close in the suggested setup, the geometry of the tissue is not important for the experiments. Therefore, simpler geometries such as a rectangular probe can be used to reduce the experimental effort. It may be possible to perform those static experiments in a petri dish instead of using a bioreactor. In doing so, multiple experiments could be executed in parallel. Using the simplified experimental setup of a rectangular probe in a static nutrition solution, correlation analyses can be performed as well. While in this work

119

the model parameters were assumed to be uncorrelated, correlations cannot be ruled out. For example, in a tissue with increased growth rates of the components, the components may decay faster as well. Furthermore, correlations between the growth rates of collagen and elastin could exist, because cells, whose productivity is above average, might synthesize both proteins faster. 8. CONCLUSION The simulation of the presented model indicates that the model is a suitable first attempt for an initial parameter analysis. Based on the performed sensitivity analysis of the model, requirements on experiments can be determined: Corresponding experiments should analyze all model parameters, but they should have special focus on the diffusion coefficient DV and the diffusion parameter γ. Future research will therefore focus on the model-based design and the execution of corresponding experiments. In turn, the model can then be evaluated and improved. All in all, a suitable initial model of the maturation process, which is easily expandable and adjustable, was developed and a sensitivity analysis method, which can not only be applied to this initial model, but can also be used to analyze prospective improved models, was performed. REFERENCES Abel, D., Bollig, A., 2006. Rapid Control Prototyping. Springer-Verlag GmbH. Albert, P. J., Schwarz, U. S., Jun. 2014. Dynamics of cell shape and forces on micropatterned substrates predicted by a cellular potts model. Biophysical Journal 106 (11), 2340–2352. Albert, P. J., Schwarz, U. S., Apr. 2016. Dynamics of cell ensembles on adhesive micropatterns: bridging the gap between single cell spreading and collective cell migration. PLoS Computational Biology 12 (4). Alberts, B., Johnson, A., Lewis, J., Morgan, D., Raff, M., Roberts, K., Walter, P., 2015. Molecular biology of the cell. Garland Science Publishing. Angello, J. C., Pendergrass, W. R., Norwood, T. H., Prothero, J., Jul. 1987. Proliferative potential of human fibroblasts: An inverse dependence on cell size. Journal of Cellular Physiology 132 (1), 125–130. Bard, J. B., Hay, E. D., Nov. 1975. The behavior of fibroblasts from the developing avian cornea. morphology and movement in situ and in vitro. The Journal of cell biology 67, 400–418. Bashey, R. I., Torii, S., Angrist, A., 1967. Age-related collagen and elastin content of human heart valves. Journal of Gerontology 22 (2), 203–208. Benton, J. A., Fairbanks, B. D., Anseth, K. S., dec 2009. Characterization of valvular interstitial cell function in three dimensional matrix metalloproteinase degradable PEG hydrogels. Biomaterials 30 (34), 6593–6603. Bouten, C. V. C., Smits, A. I. P. M., Baaijens, F. P. T., Mai 2018. Can we grow valves inside the heart? perspective on material-based in situ heart valve tissue engineering. Frontiers in Cardiovascular Medicine 5. Carmeliet, P., Jain, R. K., Sep. 2000. Angiogenesis in cancer and other diseases. Nature 407 (6801), 249–257. Crank, J., Nicolson, P., Dec. 1996. A practical method for numerical evaluation of solutions of partial differential

120

Kirsten Voß et al. / IFAC PapersOnLine 52-26 (2019) 113–120

equations of the heat-conduction type. Advances in Computational Mathematics 6 (1), 207–226. Cukier, R. I., Levine, H. B., Shuler, K. E., Jan. 1978. Nonlinear sensitivity analysis of multiparameter model systems. Journal of Computational Physics 26 (1), 1–42. Fisher, R. A., Jun. 1937. The wave of advance of advantageous genes. Annals of Eugenics 7 (4), 355–369. Flanagan, T. C., Cornelissen, C., Koch, S., Tschoeke, B., Sachweh, J. S., Schmitz-Rode, T., Jockenhoevel, S., Aug. 2007. The in vitro development of autologous fibrin-based tissue-engineered heart valves through optimised dynamic conditioning. Biomaterials 28 (23), 3388–3397. Franceschini, G., Macchietto, S., oct 2008. Model-based design of experiments for parameter precision: State of the art. Chemical Engineering Science 63 (19), 4846– 4872. Henkel, T., Wilson, H., Krug, W., Dez. 2012. Global sensitivity analysis of nonlinear mathematical models – an implementation of two complementing variancebased algorithms. In: Proceedings of the 2012 Winter Simulation Conference (WSC). IEEE, pp. 154:1–154:12. Hossain, M. S., Bergstrom, D. J., Chen, X. B., 2015. Modelling and simulation of the chondrocyte cell growth, glucose consumption and lactate production within a porous tissue scaffold inside a perfusion bioreactor. Biotechnology Reports Mai, 55–62. Kolmogoroff, A., Petrovsky, I., Piscounoff, N., 1988. Study of the diffusion equation with growth of the quantity of matter and its application to a biology problem. In: Dynamics of Curved Fronts. Elsevier, pp. 105–130. MATLAB, 2018. version 9.5 (R2018b). The MathWorks Inc., Natick, Massachusetts. Moreira, R., Neusser, C., Kruse, M., Mulderrig, S., Wolf, F., Spillner, J., Schmitz-Rode, T., Jockenhoevel, S., Mela, P., 2016. Tissue-engineered fibrin-based heart valve with bio-inspired textile reinforcement. Advanced Healthcare Materials 5 (16), 2113–2121. Mulholland, L. D., Gotlieb, I. A., Apr. 1996. Cell biology of valvular interstial cells. The Canadian journal of cardiology 12, 231–236. Olsen, L., Maini, P. K., Sherratt, J. A., 1997. A mechanochemical model for normal and abnormal dermal wound repair. Nonlinear Analysis Theory, Methods & Applications 30 (6), 3333–3338. Olsen, L., Sherratt, J. A., Maini, P. K., 1995. A mechanochemical model for adult dermal wound contraction and the permanence of the contracted tissue displacement profile. Journal of Theoretical Biology 177 (2), 113–128. Pianosi, F., Sarrazin, F., Wagener, T., Aug. 2015. A matlab toolbox for global sensitivity analysis. Environmental Modelling & Software 70, 80–85. Saltelli, A., Annoni, P., Azzini, I., Campolongo, F., Ratto, M., Tarantola, S., Feb. 2010. Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index. Computer Physics Communications 181 (2), 259–270. Schoen, F. J., 2012. Mechanisms of function and disease of natural and replacement heart valves. Annual Review of Pathology: Mechanisms of Disease 7 (1), 161–183. Shakeel, M., 2013. Travelling wave solution of the fisherkolmogorov equation with non-linear diffusion. Applied

Mathematics 04 (08), 148–160. Sherratt, M. J., Jul. 2009. Tissue elasticity and the ageing elastic fibre. AGE 31 (4), 305–325. Sodhani, D., Varun Raj, R., Simon, J., Reese, S., Moreira, R., Gesch´e, V., Jockenhoevel, S., Mela, P., Stier, B., Stapleton, S. E., Aug. 2018. Artificial textile reinforced tubular aortic heart valves – multi-scale modelling and experimental validation. In: Biomedical Technology. Springer International Publishing, pp. 185–215. Zhang, K., Lu, Z., Cheng, L., Xu, F., Jul. 2015. A new framework of variance based global sensitivity analysis for models with correlated inputs. Structural Safety 55, 1–9.