In silico exploration of post-burn contraction using uncertainty quantification

Abstract

Burns can make patients’ lives quite miserable. Apart from prominent and thickened, or hypertrophic, scars, the skin may be characterized by contraction. When this contraction is so severe that the patient loses joint mobility, it is called contracture. Then a patient may have difficulties with sports or other daily activities. The consequence can be an enormous psychosocial burden for the patient. Understanding contraction mechanisms is essential to improve and optimize the treatment of contractures. This understanding can arise from clinical (in vivo) and experimental (in vitro) observations but can also be explored using mathematical models (in silico). Mathematical models describe quantitative relations and can explain specific trends and make predictions. Further, in silico models form an alternative for animal experiments. One such mathematical model is the Biomorphoelastic model for post-burn contraction [1]. This model arises from conservation laws expressed in partial differential equations on a continuous (macro) scale. We study this model’s one- and twodimensional counterparts. The biomorphoelastic model for post-burn contraction has multiple steady states, of which some are unstable for specific parameter values. Therefore, we need a control system to stabilize and control the process at and around the steady states. In this work, we predict stability around steady states using linearized models and apply multivariate analysis to assess the local behavior around these steady states. This way, we avoid the parametric dependence of stable and unstable solutions and explain the a priori behavior of the solution. The results show that so-called signaling molecules (growth factors, cytokines, chemokines) involved in wound healing should decay quickly enough for the model to be stable. Slow signaling molecule decay results in lasting signaling, causing continued differentiation of cells that contract the tissue and increased collagen deposition that we link to hypertrophic scarring. Furthermore, the model’s convergence is monotonic, except for low viscosity values, and the numerical method has convergence order at least O(h2). There is much variation in the physiological properties of patients, therefore, giving an enormous uncertainty in post-burn contraction. One wonders, for example, why one patient develops a severe contracture and another does not. For the modeling, this is expressed in patient-specific values of the model’s input parameters, which give a patientspecific set of output variables. A significant challenge is that many patient-specific input parameter values still need to be discovered and better or more consistently documented in the literature. After all, measurements always contain uncertainty or a margin of error. It is, therefore, necessary to include uncertainty in the modeling. We assume statistical distributions for the model’s input parameters and base these on data from the literature and intuitive arguments. First, we explore the uncertainty in these input parameters by performing a sensitivity analysis. With the sensitivity analysis, we explore to what extent the outcome variables are subject to change when parameter values are varied and whether the model can make patient-specific contraction predictions. The analysis shows that parameter sensitivity significantly depends on cell proliferation rates. If the proliferation rates are equal for different cells, then the equilibrium collagen concentration significantly affects the outcome variables. This result leads to the discussion of whether we need to make the distinction between different collagen types. The analysis using unequal cell proliferation rates draws attention to the actual contraction mechanism (on a cellular level), and it discusses to what extent tissue contraction cells (myofibroblasts) proliferate. Further, the analysis provides advice for therapies. One should focus on the dermal cells and collagen to limit contraction during ‘early’ healing. However, to limit scar contraction, one should focus on post-burn signaling. Further, by performing a feasibility study, we show that the model can show significantly different contraction intensities for different groups of parameter values. The groups are chosen to relate to patients of different ages. A single simulation only provides a possible insignificant scenario from a probabilistic point of view. Therefore, we perform extensive simulations to estimate the mean outcomes, spreads, and statistical distributions for courses of contraction and patient discomfort in the first year after the burn. By drawing samples from the input parameters’ statistical distributions, we perform Monte Carlo simulations to further explore these input parameters’ uncertainty. With these Monte Carlo simulations, we determine correlations between the outputs and inputs and gain further insight into the reciprocal dependencies. Through the computations, we treat the output statistically by using histograms to model probability density functions, which provide graphical insight into the statistical distribution of the output variables. We then statistically test whether the outcome follows a particular probability distribution. For example, we use the output data to estimate the probability that the post-burn contraction will be more than a certain percentage. This uncertainty quantification method is a significant amount of work translating into computer-intensive computations, and, therefore, we approach the output variables with Bayesian finite-element trained machine learning as an alternative. For this, we use relatively small neural networks that give over 99.5% goodness of fit and spectacular computational accelerations (19354X speedup in 1D and 1815000X speedup in 2D) compared to the numerical approach. We further illustrate the use of these neural networks in an ‘online’ application. This application shows how clinicians can be offered immediate and quick access to finite element simulations and quick estimation of the probability distribution of the severity of post-burn contraction. In addition to that, this application makes such estimations accessible to clinicians. In short, we use a mathematical model that simulates post-burn contraction and patient discomfort. We study the stability of this model and quantify the uncertainties in input parameter values through sensitivity analyses and Monte Carlo studies. Because Monte Carlo simulations are computationally intensive, we approximate the output variables with neural networks and demonstrate its application. The model is stable when signaling molecules disappear from the tissue quickly enough. Furthermore, it is discussed whether different collagen types should be added to the model and to what extent tissuecontracting myofibroblasts proliferate (cell division).

Brandwonden kunnen het leven van patiënten behoorlijk ellendig maken. Afgezien van prominente en verdikte of hypertrofische littekens, kan de huid worden gekenmerkt door contractie. Wanneer deze samentrekking zo hevig is dat de patiënt de gewrichtsmobiliteit verliest, wordt dit contractuur genoemd. Dan kan een patiënt moeite hebben met sporten of andere dagelijkse bezigheden. Het gevolg kan een enorme psychosociale belasting voor de patiënt zijn. Het begrijpen van samentrekkingsmechanismen is essentieel om de behandeling van contracturen te verbeteren en te optimaliseren. Dit begrip kan voortkomen uit klinische (in vivo) en experimentele (in vitro) waarnemingen, maar kan ook worden onderzocht met behulp van wiskundige modellen (in silico). Wiskundige modellen beschrijven kwantitatieve relaties en kunnen specifieke trends verklaren en voorspellingen doen. Verder vormen in silico modellen een alternatief voor dierproeven. Een zo’n wiskundig model is het Biomorpho-elastic model for post-burn contraction [1]. Dit model komt voort uit behoudswetten uitgedrukt in partiële differentiaalvergelijkingen op continue (macro)schaal. We bestuderen de een- en tweedimensionale tegenhangers van dit model. Het biomorfo-elastische model voor contractie na verbranding heeft meerdere stabiele toestanden, waarvan sommige onstabiel zijn voor specifieke parameterwaarden. Daarom hebben we een controle systeem nodig om het proces in en rond de stationaire toestanden te stabiliseren en te beheersen. In dit werk voorspellen we stabiliteit rond stabiele toestanden met behulp van gelineariseerde modellen en passen we multivariate analyse toe om het lokale gedrag rond deze stabiele toestanden te beoordelen. Op deze manier vermijden we de parametrische afhankelijkheid van stabiele en onstabiele oplossingen en verklaren we het a priori gedrag van de oplossing. De resultaten laten zien dat zogenaamde signaalmoleculen (groeifactoren, cytokines, chemokines) die betrokken zijn bij wondgenezing snel genoeg zouden moeten vervallen om het model stabiel te houden. Langzaam verval van signaalmoleculen resulteert in blijvende signalering, wat leidt tot een voortdurende differentiatie van cellen die het weefsel samentrekken en een verhoogde collageenafzetting die we in verband brengen met hypertrofische littekens. Bovendien is de convergentie van het model monotoon, behalve voor lage viscositeitswaarden, en heeft de numerieke methode een convergentievolgorde van ten minste O(h2). Er is veel variatie in de fysiologische eigenschappen van patiënten, wat een enorme onzekerheid geeft in de contractie na de brandwond. Men vraagt zich bijvoorbeeld af waarom de ene patiënt een ernstige contractuur ontwikkelt en de andere niet. Voor de modellering wordt dit uitgedrukt in patiënt-specifieke waarden van de invoerparameters van het model, die een patiënt-specifieke set uitvoervariabelen opleveren. Een belangrijke uitdaging is dat veel patiënt-specifieke invoerparameterwaarden nog moeten worden ontdekt en beter of consistenter moeten worden gedocumenteerd in de literatuur. Metingen bevatten immers altijd onzekerheid of een foutmarge. Het is daarom noodzakelijk om onzekerheid in de modellering op te nemen. We gaan uit van statistische verdelingen voor de invoerparameters van het model en baseren deze op gegevens uit de literatuur en intuïtieve argumenten. Eerst onderzoeken we de onzekerheid in deze invoerparameters door een gevoeligheidsanalyse uit te voeren. Met de sensitiviteitsanalyse onderzoeken we in hoeverre de uitkomstvariabelen aan verandering onderhevig zijn als parameterwaarden variëren, en of het model patiëntspecifieke contractievoorspellingen kan doen. De analyse laat zien dat parametergevoeligheid in belangrijke mate afhangt van cel-proliferatiesnelheden. Als de proliferatiesnelheden voor verschillende cellen gelijk zijn, heeft de evenwichtscollageenconcentratie een significante invloed op de uitkomstvariabelen. Dit resultaat leidt tot de discussie of we onderscheid moeten maken tussen verschillende soorten collageen. De analyse met behulp van ongelijke cel-proliferatiesnelheden vestigt de aandacht op het eigenlijke contractiemechanisme (op cellulair niveau) en bespreekt in welke mate weefselcontractiecellen (myofibroblasten) prolifereren. Verder geeft de analyse adviezen voor therapieën. Men moet zich concentreren op de dermale cellen en het collageen om contractie tijdens ‘vroege’ genezing te beperken. Om littekencontractie te beperken, moet men zich echter concentreren op signalering (van signaalmoleculen) na het verbranden. Verder laten we, door een haalbaarheidsstudie uit te voeren, zien dat het model significant verschillende contractie-intensiteiten kan laten zien voor verschillende groepen parameterwaarden. De groepen zijn zo gekozen dat ze betrekking hebben op patiënten van verschillende leeftijden. Een enkele simulatie geeft slechts een mogelijk onbelangrijk scenario vanuit een probabilistisch oogpunt. Daarom voeren we uitgebreide simulaties uit om de gemiddelde uitkomsten, spreidingen en statistische verdelingen voor het verloop van contractie en ongemak voor de patiënt in het eerste jaar na de brandwond te schatten. Door steekproeven te trekken uit de statistische verdelingen van de invoerparameters, voeren we Monte Carlo-simulaties uit om de onzekerheid van deze invoerparameters verder te onderzoeken. Met deze Monte Carlo-simulaties bepalen we correlaties tussen de uit- en invoer van het model en krijgen we meer inzicht in de onderlinge afhankelijkheden. Door de berekeningen behandelen we de uitvoervariabelen statistisch door histogrammen te gebruiken om kansdichtheidsfuncties te modelleren, die grafisch inzicht geven in de statistische verdeling van de uitvoervariabelen. Vervolgens testen we statistisch of de uitkomst een bepaalde kansverdeling volgt. We gebruiken de uitvoergegevens bijvoorbeeld om de kans in te schatten dat de contractie na het verbranden meer dan een bepaald percentage zal zijn. Deze methode voor het kwantificeren van onzekerheid is een aanzienlijke hoeveelheid werk die zich vertaalt in computer intensieve berekeningen, en daarom benaderen we de uitvoervariabelen met Bayesiaanse eindige-elementen getrainde machine learning als alternatief. Hiervoor gebruiken we relatief kleine neurale netwerken die meer dan 99,5% fitheid en spectaculaire rekenversnellingen opleveren (19354x versnelling in 1D en 1815000x versnelling in 2D) in vergelijking met de numerieke benadering. We illustreren verder het gebruik van deze neurale netwerken in een ‘online’ applicatie. Deze applicatie laat zien hoe clinici onmiddellijk en snel toegang kunnen krijgen tot eindige-elementensimulaties en een snelle schatting van de waarschijnlijkheidsverdeling van de ernst van huidcontractie. Bovendien maakt deze applicatie dergelijke schattingen toegankelijk voor clinici. Kortom, we gebruiken een wiskundig model dat contractie na het verbranden en ongemak voor de patiënt simuleert. We bestuderen de stabiliteit van dit model en kwantificeren de onzekerheden in invoerparameterwaarden door middel van gevoeligheidsanalyses en Monte Carlo-onderzoeken. Omdat Monte Carlo-simulaties rekenintensief zijn, benaderen we de uitvoervariabelen met neurale netwerken en demonstreren we de toepassing ervan. Het model is stabiel als signaalmoleculen snel genoeg uit het weefsel verdwijnen. Verder wordt besproken of verschillende soorten collageen aan het model moeten worden toegevoegd en in hoeverre weefsel contracterende myofibroblasten prolifereren (celdeling).