UvA-DARE (Digital Academic Repository) Inverse uncertainty quantification of a mechanical model of arterial tissue with surrogate modelling

Disorders of coronary arteries lead to severe health problems such as atherosclerosis, angina, heart attack and even death. Considering the clinical significance of coronary arteries, an efficient computational model is a vital step towards tissue engineering, enhancing the research of coronary diseases and developing medical treatment and interventional tools. In this work, we applied inverse uncertainty quantification to a microscale agent-based arterial tissue model, a component of a three-dimensional multiscale in-stent restenosis model. Inverse uncertainty quantification was performed to calibrate the arterial tissue model to achieve a mechanical response in line with tissue experimental data. Bayesian calibration with a bias term correction was applied to reduce the uncertainty of unknown polynomial coefficients of the attractive force function and achieve agreement with the mechanical behaviour of arterial tissue based on the uniaxial strain tests. Due to the high computational costs of the model, a surrogate model based on the Gaussian process was developed to ensure the feasibility of the computations.


Introduction
Cardiovascular diseases (CVDs) are the top cause of death globally, leading to an estimated 16% of total deaths annually [1].Atherosclerosis, a narrowing or blocking of arteries due to a plaque buildup is a common CVD, which is often treated by performing percutaneous coronary interventions (PCIs) with stenting in coronary vessels [2,3].However, an incidence rate of an excessive amount of neointima formation after PCI and a repeat narrowing of the vessel, a condition known as in-stent restenosis (ISR), is notably high, making it a public health concern of significant importance [4,5].A three-dimensional computational model ISR3D has been developed to simulate smooth muscle cells proliferation and the ISR process, as well as to predict the restenosis progression under different scenarios [6][7][8].
The ISR3D is a coupled multiscale model consisting of three singlescale submodels: an initial condition model, a vessel tissue model, and a blood flow model, as well as utility modules that facilitate communication between submodels [6].The vessel tissue submodel is a cell-scale agent-based model (ABM), and it uses pairwise repulsive and attractive forces between cell-representing agents to model mechanical behaviour.Three layers of arterial tissue -diseased intima with atherosclerosis, media and adventitia are modelled explicitly.The microscale interactions result in the mechanical response representative of biological tissue on a macroscale.The attractive force incorporates macroscopic properties into the cell-scale model.It is modelled as a 6th order polynomial pairwise function, depending on cells' size and bond strain.The polynomial coefficients are model-specific, and no prior values of these coefficients are known.A goal of this work is to calibrate those unknown polynomial coefficients to match the macroscopic behaviour (i.e. the stress-strain experimental data of human coronary arteries) of the arterial tissue according to [9,10].The considered macroscopic data is a stress-strain relation for uniaxial strain tests of tissue stretched in a circumferential direction, which is a simple but effective representation of the strain occurring in the tissue during stenting or pressurisation.The existing continuum models with finite element (FE) methods have managed to capture the mechanical response of arterial tissues well [11,12].Therefore, we utilised a continuum model to generate the data and subsequently applied them to the ABM calibration.By integrating continuum modelling, which provides a macroscopic view of the system with agent-based modelling, that allows for a more detailed understanding of the microscopic behaviour of individual agents on the cellular level [7,13,14], researchers can gain a better understanding of complex systems at multiple scales [15].https://doi.org/10.1016/j.ress.2023 Inverse uncertainty quantification (IUQ) is a process that determines uncertain inputs based on experimental data or a measured output using Bayesian techniques [16,17].It allows people to integrate prior information of the uncertain inputs and update that knowledge with existing data based on the Bayes' rule [18,19].Since the marginal likelihood, used for normalising the posterior distribution, is typically computationally intractable, multiple approximation methods have been developed to estimate the posterior distributions or the samples of the posterior distribution, such as Markov chain Monte Carlo (MCMC) [20,21] and variation inference [22].The derived posterior distribution can be subsequently applied for predictive distribution.As a result, IUQ has found a wide application in the field of computational science and engineering, particularly for solving inverse problems [23][24][25][26].
IUQ typically requires a large number of model evaluations.Such many-query scenarios may become computationally prohibitive for a computationally expensive model.The problem can be addressed by employing surrogate modelling, that approximates the original complex model.One of the main categories of surrogates is a group of data-driven methods, which considers the model as a black box and empirically approximates the latent mapping between inputs and model responses.Typical instances of data-driven methods include radial basis function method [27], polynomial chaos expansion [28][29][30], neural networks [31,32], support vector machines [33], et cetera.The Gaussian process is another state-of-the-art regression method and is widely applied due to its non-parametric and probabilistic nature [34][35][36].It assumes that the responses of the model follow multivariate normal distributions and, conditioning on existing data, a response at a new input position can be predicted.Its hyperparameters can be learned through maximum marginal likelihood.We developed the Gaussian process regression surrogate to substitute the original ABM for the evaluations in the inverse uncertainty quantification process.
In this paper, we developed a Bayesian framework combining inverse uncertainty quantification techniques and a Gaussian process surrogate model to calibrate a microscale agent-based model of arterial tissue based on a macroscale continuum constitutive model for the mechanical behaviour of arterial tissue.This approach improved the model's accuracy, which serves as a submodel in the multiscale simulation of a complex process of in-stent restenosis.The reliability and accuracy of each component have the utmost importance for the overall reliability and precision of the complex system driven by the integration and interaction of the component modules.
The paper is arranged as follows.The material, agent-based and surrogate models of arterial tissue are described in Section 2. In the same section, the IUQ process is as well introduced.The result of surrogate modelling and IUQ are presented in Section 3, followed by the discussion and conclusion in Sections 4 and 5.The nomenclature used in the paper is given in Table 1.

Agent-based model of arterial tissue
The arterial tissue model is a part of the 3D multiscale model of in-stent restenosis (ISR3D), extensively described in [7].From the mechanical point of view, the vessel wall is modelled as a centre-based agent-based model, where individual cells are modelled as spheres, interacting via adhesive and repulsive forces.Spheres have their radii set to match the volume of the cells, and forces act between cell centres.In ISR3D, biological behaviour is modelled by imposing a biological ruleset on each agent, while simplified mechanical properties of the centre-based approach allow massive simulations.The ABM of the arterial wall consists of diseased intima with atherosclerosis, media and adventitia.Fig. 1(a) illustrates the three-layered agent-based vessel wall with a placed stent.For each layer of the ABM of arterial tissue, isotropic tissue with an almost constant density is assumed, which is generated by placing cells randomly, maintaining a minimum distance between neighbouring pairs using a three-dimensional Poisson disc sampling by Bridson's algorithm [37].Interactions occur between adjacent cells, and neo-Hookean and polynomial attractive forces are selected for repulsion and adhesion, respectively [7].
It should be noted that an actual arterial tissue is highly anisotropic.However, we choose to focus on an isotropic formulation as our primary goal is to apply the tissue model to stenting modelling and vessels reconstructed from, e.g.optical coherence tomography or intravascular ultrasound images [38].Thus, most of the stress in our use scenarios will likely be in the circumferential direction.Also, finding the fibre directions in an atherosclerotic plaque is technically challenging and not something routinely done -e.g.Akyildiz et al. [39] get a detailed image of fibres by studying a plaque ex vivo.Since this type of data is unexpected for most of the stented vessels, we opt for the model calibration based on the circumferential behaviour of arterial tissue.

Mechanical properties of tissue ex vivo
One of the most widespread techniques to determine the mechanical properties of arterial tissue is ex vivo uniaxial tensile testing.Strips of dissected tissue are subjected to uniaxial tension tests by applying an increasing force to the specimen, and the deformation and ultimate tensile stresses are measured.This provides stress-strain relation for each tested direction of the specimen.The (engineering) strain is defined as the elongation divided by the initial length, and the stress is calculated as the force divided by the cross-sectional area of the specimen [40].In this paper, engineering stress is calculated to match the experimental approach from Holzaphel et al. [10] and other papers.This way, the initial cross-sectional area is used in the stress calculations.The stress-strain relations for uniaxial strain tests of tissue stretched in a circumferential direction effectively represent  the strain the tissue undergoes during the stenting or pressurisation procedures.The computational tissue model calibrated based on the stress-strain behaviour achieving a mechanical response representative of experimental data, can mimic the realistic macroscopic behaviour of arterial tissue.

Material model of arterial wall
The constitutive model for hyperelastic, incompressible rubber-like material is used to produce ground-truth stress-strain relations of uniaxial tension tests of the arterial tissue for the further calibration task.The generalised polynomial form of the strain energy function is used, which was introduced in the seminal work of Rivlin and Saunders [41]: where  is a strain energy function and describes a material response while   is a material parameter that explains the material's shear behaviour.
are two strain invariants, and   ,  = 1, … , 3 is a principal extension ratio, which is in the following relation with a principal strain value   :   = 1 +   .The 6th order reduced polynomial strain energy function was employed, which is derived by setting  = 0 in Eq. ( 1), resulting in the reduced strain energy function to the first strain invariant.The rationale behind reducing the general polynomial function is reviewed by Lapeer et al. [42].The reduced form of the general polynomial strain energy function has a form [43]: where  = 6 is the order of the reduced polynomial energy function.
Since we aim to obtain the stress-strain relation for the uniaxial strain deformation that acts in a single direction, the following associations hold for the principal stretches and stretch and strain in the loading direction - and  [41]: The nominal stress for the uniaxial tensile (UT) deformation is obtained by the formula [44]: As a result, the nominal stress-strain relationship based on the 6th order reduced polynomial strain energy function can be derived: Fig. 2 shows stress-strain curves obtained by the material model for three layers of arterial tissue.The choices of the 6th order reduced polynomial strain energy function and the material parameters  0 , demonstrated in Table 2, are according to the FE arterial model by Zun et al. [9], which achieved a macroscopic behaviour of an isotropic macroscale model approximating the range of variability of experimental data from a human coronary artery by Holzaphel et al. [10].Thus, the stress response of the material model, later used for the model calibration in the IUQ process, is an analytical isotropic approximation of the experimental in vivo data.

Mechanical properties in agent-based arterial tissue model
Attractive forces incorporate multiple methods of cell-cell interaction into the model, such as extracellular fibres or cell adhesion.6th order polynomial pairwise interaction force is used as an agent-agent attraction bond force [9]: 1 and  2 are the radii of two interacting cells,  is a bond strain and   , for  = 1, … , 6 are coefficients of the attraction bond force.
The polynomial coefficients are model-specific and do not have any prior meaning or values.The aim is to find those coefficients for the attraction bond force to achieve the macroscopic behaviour of arterial tissue.The bond strain  is quantified by the formula: Where  is the distance between the centres of two interacting cells.The choice of the 6th order follows the order of the analytical reduced polynomial strain energy function.The corresponding force coefficients are constrained to be non-negative to achieve a minimum of potential energy for each bond at  = 0.For different layers of arterial tissue, separate sets of force coefficients and, thus, distinct bond forces are used.Uniaxial strain tests in the agent-based arterial tissue model are performed separately on intima, media and adventitia in the following manner: a single side of the tissue is fixed, while the rest is stretched up to a particular strain value, and the opposite side is set as well.The system of agents is then solved to find a final shape of the tissue with these boundary conditions.Afterwards, the forces applied to the fixed cells are summed up and divided by the cross-section area to measure the stress for the strain value.Fig. 1(b) presents a generated tissue sample of a single layer in an unstrained state, while Fig. 1(c) shows a tissue that is subjected to the uniaxial strain, with force coefficients  1 = 2.14 and  2 , … ,  6 = 0.The values are selected to obtain a simple linear fit to the high-strain part of the media data shown in Fig. 2. The microscopic behaviour of the arterial tissue model strongly depends on the tissue structure, such as the radii of cells, their arrangement, and density of the tissue.A detailed investigation of the effect of the density of the tissue on its macroscopic behaviour is presented in the supplementary material -Appendix A.

Surrogate model
Here, a Gaussian process regression surrogate model is introduced to represent the latent function   that maps input strain values and force coefficients to output stress values.Assume that a stress value  ∈ R can be considered as a function of a strain value and polynomial coefficients of attractive force  = ( 1 , …  6 , ) ∈ R 7 , for: is a Gaussian noise term, describing stochasticity in observations with a variance of  2  .Given a set of evaluated inputs and outputs of size   , {(, )} = {(  ,   )}   =1 , the covariance of  can be written as: where (, ) is a covariance matrix describing correlations between function values at different input points.A Gaussian process is a collection of random variables, any finite collection of which follows a joint multivariate normal distribution [45]: Where  and  denote the mean and the covariance functions, respectively, the mean function is generally set to zero to avoid expensive computations.The covariance function, also known as a kernel, contains hyperparameters such as length-scale, signal variance and noise variance.They regulate a priori correlation between arguments.Based on the kernel selection procedure presented in the supplementary material -Appendix B -the Matern kernel is selected as a covariance function [46].The hyperparameters are fine-tuned by optimising the log marginal likelihood function, which has the following form [45]: Subsequently, evaluations of new unobserved data points are predicted from the resulting posterior distribution.Initially, prior hyperparameters of the GP model kernel are set to 1 and are optimised using the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm [47].The GP regression model is implemented via Gaussian Process framework GPy [48].

Performance evaluations
To evaluate the predictive capability of a surrogate model, we assess the prediction accuracy of the surrogate via the root mean square error (RMSE), standardised root mean squared error (SRMSE) and the coefficient of determination (also referred to as  2 ). 2 is a ratio of the output variance explained by the surrogate model [49], and its value ranges from 0 (null predictivity) to 1 (perfect predictivity).

Experimental design
Separate layers of arterial tissue have different mechanical properties [10].They are represented by three distinct interaction forces and are calibrated individually.Thus, separate sets of coefficients are found for intima, media and adventitia.
ABM stress responses for uniaxial strain tests on an extensive collection of sets of coefficients were obtained to gather sufficient data for building the surrogate model.A uniaxial strain test is performed by stretching the tissue to the strain value varied from 0 to 0.35 with an increment of 0.005, giving a total number of observations   = 71 of stress-strain sets.
Table 3 shows preliminary sets of force coefficients as an initial guess for intima, media and adventitia layers, such that when given to the ABM, they produce stress-strain curves in a rough agreement with the benchmark data in terms of magnitude.The ranges of coefficients for three layers and ''Joint'' bounds covering all three ranges are demonstrated in Table 3.All the upper and lower bounds for each layer are selected around the preliminary guess, which is justified by the sensitivity analysis of each coefficient, demonstrated in the supplementary material -Appendix C.
Sets of coefficients  = ( 1 , …  6 ) were sampled from the designed ranges using Latin Hypercube Sampling (LHS).LHS, used to achieve good parameter space coverage and avoid clustering, was implemented via an open-source Python package SMT [50].The validation data contains a large number of input parameter sets for all strain values, resulting in stress outcomes of different magnitudes, which would result in high sensitivity of error measures to the scale of the target values.To tackle this issue, the RMSE was calculated for each strain value separately and was normalised by the variance of the target stress values [45], giving the SRMSE values, which were finally summed up to obtain the total SRMSE term for each layer.

Inverse uncertainty quantification 2.6.1. Bayesian calibration
Bayesian calibration relates prior information with uncertainty to posterior information based on the likelihood of simulated outputs from the computational model.In each iteration of the Bayesian calibration, the posterior probability density functions of calibration parameters are updated in a way that is most likely to align with the benchmark data.The posterior distributions of the parameters inferred after a sufficient number of iterations are the most probable calibration of the model, which means that a model calibrated with Bayesian calibration can efficiently produce expected behaviour by determining the best estimate of parameter uncertainties.
We aim to calibrate the set of attractive force parameters  based on the evidence data of stress output  during the stress-strain tests.Bayesian calibration is an application of the Bayesian inference method, in which the probability for a hypothesis is updated as more evidence is provided: where () is a prior distribution of the model parameters , (|) is known as likelihood or sampling distribution and () = ∫  ()(|) is integrated over the full parameter space .However, since () does not depend on  and for fixed  it can be considered as a constant, equation ( 12) can be reformulated as unnormalised posterior density [51]:

IUQ model formulation
The iterative Bayesian calibration process is conducted based on the model updating equation -''true value = simulated value + uncertainty'', which links observation values to the model response and uncertainty and is used to tune unknown model calibration parameters and the discrepancy term.In this work, a model updating equation without a model discrepancy function was selected since combining calibration of the model parameters and the discrepancy function faces the non-identifiability problem due to the challenging task of jointly fitting the model and model discrepancy.The term refers to the state when the model choice and the selected calibration targets are insufficient to obtain the unique values of the calibration parameters [52].Besides, in most cases, it is difficult to distinguish the effects of the discrepancy function and calibration parameters on response predictions when both are incorporated in the calibration process [53].We introduce a zeromean Gaussian bias term that accounts for aleatory and epistemic uncertainty sources, and its variance is calibrated throughout the IUQ process [54,55].
Inverse uncertainty quantification is formulated using the terminology of Kennedy and O'Hagan [56].The model has variable and calibration inputs.The variable inputs have known values for each observation used in IUQ, while the calibration inputs are the unknown parameters we aim to identify.The variable inputs  = ( 1 , … ,    ) are   observed strain values, obtained by stretching the tissue so that the corresponding strain value varies from 0 to 0.35 by an increment of 0.025.  = ( 1, , … ,  6, ),  = 1, … ,   is a vector of unknown attractive force coefficients we aim to calibrate, where   is a size of sampled coefficient space, consisting of sets of different combinations of coefficient values.
The benchmark data used for calibration is the analytical stress response of uniaxial stress-strain tests, obtained by the material model and is denoted by  = ( 1 , … ,    ).The IUQ model stress responses for sets of the input parameters (  , ), for  = 1, … ,   are produced as predictions of the Gaussian process surrogate model, denoted by  = ( 1 , … ,    ), where   is a total number of predictions.The set of stress data for the IUQ process is then denoted by  = (, ).As a result, the following relation is inferred between the analytical values of stress and the ABM response: where   -describes the original arterial tissue ABM we are calibrating,  is a total uncertainty of the IUQ model, while  * is the reference value of calibrated parameters and the quantity of interest, such that given to the ABM, stress response for the uniaxial strain tests aligns with the analytical solution.

Modelling prediction uncertainty
The  term represents a total prediction uncertainty, which can be categorised into the four main groups: observation error, residual variability, code uncertainty and model inadequacy [56].The observation error is the uncertainty caused by the measurement noise when collecting the experimental data.At the same time, the residual variability is the aleatory uncertainty of the model when parameters and conditions are fully specified and fixed.The code uncertainty comes from the fact that ABM response at any given set of inputs used in IUQ cannot be obtained due to the high computational costs; thus, metamodel for predicting stress value in unknown input points is introduced, which raises code or also known as  .The model inadequacy also called   can be caused by missing physics, inaccurate modelling assumptions, numerical errors or any other causes that cannot be assessed in advance in contrast with the error terms mentioned earlier.Since our calibration data is a simulation output of the material model of arterial tissue and not measurements from in vivo or in vitro experiments, observation error can be discarded.We only consider residual variability   , the model inadequacy   and the code uncertainty   .
The total uncertainty is broken into three uncertainty components and is represented as their sum, following De Vries et al. [55]: Uncertainty terms are assumed to follow independent Gaussian distributions with zero mean, and corresponding covariance matrices: Besides, the error variables are assumed to be independently and identically distributed.Thus, they are uncorrelated and their covariance matrices have diagonal structures, with the diagonal elements representing the variances and the off-diagonal elements equal to zero.Due to the independent Gaussian distribution of uncertainties, the total prediction error term  has the Gaussian distribution and its covariance matrix has the following form: As mentioned above,   and   are known a priori.Since stochasticity of the ABM or GP prediction uncertainty is not prevalent, for simplicity, homoscedastic error terms with equal variances are assumed: where  is an identity matrix,  2  is a mean variance of 100 simulations of the ABM stress response for uniaxial strain tests per strain values, while  2   is an averaged GP's prediction variance for the validation data.The variance parameter  2  , used for modelling uncertainties in   , is estimated along the IUQ process.As a result, the full set of uncertain calibration parameters is defined as  = (,  2  ), where  and  2  are assumed to be independently distributed.The likelihood is modelled as a joint Gaussian distribution [45].Its mean is the calibration data given the GP model response for the prior uncertain parameters, and the covariance matrix   is from the total uncertainty : ) Since optimal parameter ranges were not known in advance, uniform distribution was used to obtain the prior marginal distributions.The lower and upper bounds of uniform distribution were according to Table 3.

Calibration process
After defining prior distributions for the calibration parameters and the likelihood function, a MCMC sampling with a Metropolis-Hastings algorithm was used to sample from the obtained posterior distributions [57].An initial state for the sampler in the parameter space by the maximum a posteriori (MAP) method was found, which is a numerical optimisation method to find a point estimate of the mode of the distribution [58].The calibration process was then performed iteratively to update our beliefs about the calibration parameters.The obtained posterior distributions were used as the following inference's prior distributions, and new posteriors were produced.This process of the calibration was then iterated until the posterior distributions converged.For reusing posteriors, first, Gaussian kernel density estimation (KDE) [59] was used to estimate the PDF of a random variable in a non-parametric way.To ensure convergence to true parameter values, independent data is needed in each iteration.Thus, for sampling from obtained PDFs, a linear interpolation of PDF was evaluated on evenly distributed points on the extended domain of the posterior samples.The Bayesian calibration process was performed using a probabilistic programming package for Python -PyMC3 [60].

Performance of GP
A single Gaussian process model was trained on the parameter space constructed with the   = 15 strain values and 1000 sets of force coefficients sampled from the ''Joint'' ranges given in Table 3.The ratio between a training and validation set was 70% to 30%.Gaussian process predictions were compared to the ABM stress responses on the validation data for the complete set of strain values consisting of   = 71 observations; thus, interpolation performance in unknown strain values was also monitored.The quality of the surrogate predictor was assessed via the total SRMSE and the coefficient of determination  2 on the validation data.As demonstrated in Table 4, the values of total SRMSEs, calculated per sets of coefficients from individual sampling ranges of each layer, are low, and  2 coefficients are close to 1.0, which suggests satisfactory predictive capabilities of the surrogate.The investigation of the mean difference between GP prediction with a 95% confidence interval (CI) and ABM stress response is presented in the supplementary material -Appendix D.
For an intelligible illustration of the predictive capabilities of the surrogate, stress-strain curves produced by the ABM and GP model predictions given the set of preliminary coefficients and the full set of strain values were compared to each other.Fig. 3 illustrates this comparison for intima, media and adventitia layers.In all three cases, the prediction of the surrogate model is in line with ABM stress behaviour, producing low RMSEs and narrow confidence intervals of predictions.

Inverse uncertainty quantification
Inverse uncertainty quantification, using Bayesian calibration of the polynomial force coefficients and the model uncertainty term, was conducted for three layers of the arterial tissue.Fig. 4 illustrates PDFs of model parameters during the Bayesian calibration process.PDFs from earlier calibration steps are nearly flat, close to a prior uniform distribution, referring to flat marginal likelihoods of the model parameters, which means that there is a big number of maximum likelihood estimates of parameters.In the later iterations, we see a reduction in uncertainty and posterior PDFs that turn into narrow probability distributions, far from the prior, signifying model identifiability.Finally, a convergence of marginal posterior distributions of the parameters to the targeted distributions becomes apparent for three layers after 80 iterations of Bayesian calibration.The calibration process on the plots is presented via the sequential colour scheme, where initial PDFs are shown by the shades of yellow colour, which are monotonically getting darker throughout the calibration process and finally converge to targeted posterior distributions presented with the dark blue colour.
Table 5 shows the statistics of the converged posterior distributions of 6 attractive force coefficients and the variance of the model inadequacy term after 80 iterations for three layers of the tissue.More  specifically, the mean and standard deviation of the resulted probability distributions.Besides, the potential scale reduction factors (PSRF) were calculated, which can be considered as a convergence diagnostic [61].PSRF values were very close to 1 for each layer and calibration parameter, indicating that associated chains likely converged to the targeted posterior distributions [62].Thus, running simulations any longer was considered unnecessary.
The marginal distribution of  2  , the variance of the model uncertainty term   , is centred around the mean of 0.001, which validates the model and the assumption that the prediction error is normally distributed with a zero mean.For assessing the precision of the IUQ results, mean values of the final marginal posterior distributions of the coefficients were used for collecting the GP predictions and ABM responses per the specified set of force coefficients and the complete set of observed strain values.Fig. 5 illustrates the obtained stress-strain curves of ABM compared to the analytical data of the material model.The ABM of arterial tissue shows sufficient precision in fitting the ground-truth uniaxial strain test data with RMSE value up to 2.1%, which verifies that IUQ procedure was successful and the aim of the ABM to replicate realistic mechanical properties was achieved.

Speed up
The ABM of arterial tissue is implemented using OpenMP, which supports shared-memory, parallel multiprocessing programming.The uniaxial strain tests were run on the 16-core node of the SURFsara Lisa cluster in parallel, utilising all 16 cores for each test.The average simulation time of a single uniaxial strain test ranges between 1.3 to 1.5 min wall clock and depends on the scale of the force coefficient values, taking longer time for the higher values and vice versa.One Bayesian calibration iteration of the IUQ for satisfactory outcomes requires at least 4000 draws of the model stress response.Thus, a single calibration iteration for one out of three arterial layers would take 96 h if running the original model.On the other hand, the GP model takes 0.2 ms on average to produce the stress prediction, irrespective of the magnitude of input coefficient values.5, versus analytical data for three layers of arterial tissue.

Discussion
The benchmark data every IUQ and calibration problem relies on was obtained analytically by following an isotropic material model of arterial wall [9], which is a representation of in vitro behaviour of arterial tissue, approximating the range of variability of the experimental stress-strain curves [10].Flexibility to generate desirable data allowed for efficient calibration and validation processes.GP regression approximated the original model's stress behaviour in response to the uniaxial tissue stretching in a circumferential direction.Obtaining a precise and computationally efficient representation of the ABM's mechanical behaviour was vital for effective sensitivity analysis and calibration processes.
A Bayesian calibration with a bias term correction was performed as a technique of IUQ to quantify uncertainties of attractive force coefficients and estimate their values based on the analytical data.Since our objective was to find the precise values of the polynomial force function coefficients, the challenge of distinguishing between the effects of calibration parameters and other factors could have been limited.Besides, as the coefficients do not have biological value, the limitation of possibly converging to ''pseudo-true'' values, lacking physical interpretation, was not much of a concern.
When attractive forces were formed based on the calibrated coefficients, performing uniaxial strain tests on the ABM showed a stressstrain relation in line with analytical data for all three layers.As a result, uncertainties about the unknown parameters were reduced, and the attractive forces of the ABM of arterial tissue layers were determined, providing a reasonable macroscopic stress-strain relationship for uniaxial strain tests.The efficient computational model of arterial tissue can ensure the reliable and accurate behaviour of a larger-scale model with a broader scope and diverse applications, of which it is a component.Specifically, the macroscale mechanical model can be used to provide mechanical information to a biological model of cells on a microscale.The direct application of the obtained arterial tissue model with realistic macroscopic behaviour will simulate the implantation of a stent for the ISR3D model.The paper by Zun et al. [7] shows the application of ISR3D by modelling the In-stent restenosis (ISR) process in porcine coronary arteries and validating the results by in vivo data.Furthermore, the arterial tissue model can contribute to advancing vascular medicine and clinical treatment of artery diseases, designing interventional devices such as vascular implants [63] and limiting unfeasible, expensive and sometimes unethical in vivo or in vitro studies.A viable microscale agent-based model of arterial tissue can be adapted to other mechanobiological applications, such as models of valves or vein grafts.
A limitation of the arterial tissue model is that only the isotropic formulation is considered and it is calibrated to match the circumferential behaviour of the vessel.However, extending the model to include anisotropy is in our future plans.It should also be noted that the approach outlined in this paper can be generalised to other pairwise interaction forces, including anisotropic ones if uniaxial strain tests in the axial and radial directions are also included.The methodology presented here can be applied to other microscale models where mechanical properties are important.Furthermore, according to the complexity of the phenomenon and dynamics of interest, more elaborate methodological approaches can be incorporated, such as modelling inadequacy function when the IUQ parameters have physical values [64]; using the active learning strategy for the surrogate modelling [65]; reducing dimensionality incorporating principal components that account for uncertainty of the high-dimensional dynamic output [66]; and/or calibrating the parameters of the surrogate model in the IUQ process [67] if the pre-trained surrogate model does not achieve the sufficient precision.

Conclusions
This paper applied inverse uncertainty quantification to determine interaction forces in the cell-resolved agent-based arterial tissue model from the analytical data of biological tissue's macroscopic behaviour.Considering the computational intensity of the model, the necessity of surrogate modelling was anticipated.Overall, attractive force coefficients were successfully calibrated using the proposed IUQ model, which means that interaction forces between the agent-based arterial tissue model cells were found, providing a reasonable macroscopic stress-strain relationship for uniaxial strain tests.As a result, the model reveals a realistic mechanical behaviour of biological tissue.
The result is a versatile and generalisable approach for modelling and calibrating microscale model of natural phenomena using inverse uncertainty quantification techniques and surrogate modelling based on a macroscale mechanical model.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.(a) An illustration of the agent-based model of arterial vessel wall consisting of three layers: diseased intima (light blue), media (dark blue), adventitia (red) and a metal stent.The single layer of arterial tissue before (b) and after (c) applying a particular value of uniaxial strain with the force coefficients  1 = 2.14 and  2 , … ,  6 = 0.

Fig. 2 .
Fig. 2. Stress-strain curves obtained by the material model with the 6th order reduced polynomial form of the strain energy function for three layers of arterial tissue.

Fig. 4 .
Fig. 4. 1D marginal posterior distributions of the calibration parameters for the intima, media and adventitia layers of arterial tissue after 80 iterations of Bayesian calibration.

Fig. 5 .
Fig. 5. ABM stress response for uniaxial strain tests with the calibrated set of force coefficients, given in Table5, versus analytical data for three layers of arterial tissue.

Table 2
[9]erial model parameters according to[9], describing three layers of the arterial wall: intima, media and adventitia.

Table 3
The preliminary sets (PS) of the attractive force coefficients and their lower (LB) and upper (UB) bounds of the sampling ranges for three layers (intima, media and adventitia) of the arterial tissue model and ''Joint'' bounds covering all three ranges.

Table 4
Total SRMSE and predictivity coefficient  2 of GP surrogate model for three layers of the arterial tissue model.

Table 5
Statistics of converged posterior distributions for three layers of arterial tissue, namely, the mean and standard deviation of the resulted probability distributions after 80 iterations of Bayesian calibration.GP prediction with 95% CI versus ABM stress response with the preliminary set of force coefficients for intima, media and adventitia.