In this book we discuss several numerical methods for solving ordinary differential equations. We emphasize the aspects that play an important role in practical problems. We confine ourselves to ordinary differential equations with the exception of the last chapter in which we discuss the heat equation, a parabolic partial differential equation. The techniques discussed in the intro-ductory chapters, for instance interpolation, numerical quadrature and the solution to nonlinear equations, may also be used outside the context of differential equations. They have been in-cluded to make the book self-contained as far as the numerical aspects are concerned. Chapters, sections and exercises marked with a * are not part of the Delft Institutional Package. The numerical examples in this book were implemented in Matlab, but also Python or any other programming language could be used. A list of references to background knowledge and related literature can be found at the end of this book. Extra information about this course can be found at http://NMODE.ewi.tudelft.nl, among which old exams, answers to the exercises, and a link to an online education platform. We thank Matthias Moller for his thorough reading of the draft of this book and his helpful suggestions.@en

In this paper a Spatial Markov Chain Cellular Automata model for the spread of viruses is proposed. The model is based on a graph with connected nodes, where the nodes represent individuals and the connections between the nodes denote the relations between humans. In this way, a graph is connected where the probability of infectious spread from person to person is determined by the intensity of interpersonal contact. Infectious transfer is determined by chance. The model is extended to incorporate various lockdown scenarios. Simulations with different lockdowns are provided. In addition, under logistic regression, the probability of death as a function of age and gender is estimated, as well as the duration of the disease given that an individual dies from it. The estimations have been done based on actual data of RIVM (from the Netherlands).

@en

Severe burn injuries often lead to skin contraction, leading to stresses in and around the damaged skin region. If this contraction leads to impaired joint mobility, one speaks of contracture. To optimize treatment, a mathematical model, that is based on finite element methods, is developed. Since the finite element-based simulation of skin contraction can be expensive from a computational point of view, we use machine learning to replace these simulations such that we have a cheap alternative. The current study deals with a feed-forward neural network that we trained with 2D finite element simulations based on morphoelasticity. We focus on the evolution of the scar shape, wound area, and total strain energy, a measure of discomfort, over time. The results show average goodness of fit (R²) of 0.9979 and a tremendous speedup of 1815000X. Further, we illustrate the applicability of the neural network in an online medical app that takes the patient's age into account.@en

We consider the stability analysis of a two-dimensional model for post-burn contraction. The model is based on morphoelasticity for permanent deformations and combined with a chemical-biological model that incorporates cellular densities, collagen density, and the concentration of chemoattractants. We formulate stability conditions depending on the decay rate of signaling molecules for both the continuous partial differential equations-based problem and the (semi-)discrete representation. We analyze the difference and convergence between the resulting spatial eigenvalues from the continuous and semi-discrete problems.

@en

Cancer cell migration between different body parts is the driving force behind cancer metastasis, which causes mortality of patients. Migration of cancer cells often proceeds by penetration through narrow cavities in possibly stiff tissues. In our previous work [12], a model for the evolution of cell geometry is developed, and in the current study we use this model to investigate whether followers among (cancer) cells benefit from leading (cancer) cells during transmigration through micro-channels and cavities. Using Wilcoxon’s signed-rank text on the data collected from Monte Carlo simulations, we conclude that the transmigration time for the stalk cell is significantly smaller than for the leading cell with a p-value less than 0.0001, for the modelling set-up that we have used in this study.

@en

Deep dermal wounds induce skin contraction as a result of the traction forcing exerted by (myo)fibroblasts on their immediate environment. These (myo)fibroblasts are skin cells that are responsible for the regeneration of collagen that is necessary for the integrity of skin We consider several mathematical issues regarding models that simulate traction forces exerted by (myo)fibroblasts. Since the size of cells (e.g. (myo)fibroblasts) is much smaller than the size of the domain of computation, one often considers point forces, modelled by Dirac Delta distributions on boundary segments of cells to simulate the traction forces exerted by the skin cells. In the current paper, we treat the forces that are directed normal to the cell boundary and toward the cell centre. Since it can be shown that there exists no smooth solution, at least not in H¹ for solutions to the governing momentum balance equation, we analyse the convergence and quality of approximation. Furthermore, the expected finite element problems that we get necessitate to scrutinize alternative model formulations, such as the use of smoothed Dirac Delta distributions, or the so-called smoothed particle approach as well as the so-called ‘hole’ approach where cellular forces are modelled through the use of (natural) boundary conditions. In this paper, we investigate and attempt to quantify the conditions for consistency between the various approaches. This has resulted into error analyses in the L²-norm of the numerical solution based on Galerkin principles that entail Lagrangian basis functions. The paper also addresses well-posedness in terms of existence and uniqueness. The current analysis has been performed for the linear steady-state (hence neglecting inertia and damping) momentum equations under the assumption of Hooke's law.

@en

Health care is undergoing a profound technological and digital transformation and has become increasingly complex. It is important for burns professionals and researchers to adapt to these developments which may require new ways of thinking and subsequent new strategies. As Einstein has put it: "We must learn to see the world anew." The relatively new scientific discipline "Complexity science" can give more direction to this and is the metaphorical open door that should not go unnoticed in view of the burn care of the future. Complexity science studies "why the whole is more than the sum of the parts." It studies how multiple separate components interact with each other and their environment and how these interactions lead to "behavior of the system." Biological systems are always part of smaller and larger systems and exhibit the behavior of adaptivity, hence the name complex adaptive systems. From the perspective of complexity science, a severe burn injury is an extreme disruption of the "human body system." But this disruption also applies to the systems at the organ and cellular levels. All these systems follow the principles of complex systems. Awareness of the scaling process at multilevel helps to understand and manage the complex situation when dealing with severe burn cases. This paper aims to create awareness of the concept of complexity and to demonstrate the value and possibilities of complexity science methods and tools for the future of burn care through examples from preclinical, clinical, and organizational perspectives in burn care.

@en

Skin contraction is an important biophysical process that takes place during and after recovery of deep tissue injury. This process is mainly caused by fibroblasts (skin cells) and myofibroblasts (differentiated fibroblasts which exert larger pulling forces and produce larger amounts of collagen) that both exert pulling forces on the surrounding extracellular matrix (ECM). Modelling is done in multiple scales: agent-based modelling on the microscale and continuum-based modelling on the macroscale. In this manuscript we present some results from our study of the connection between these scales. For the one-dimensional case, we managed to rigorously establish the link between the two modelling approaches for both closed-form solutions and finite-element approximations. For the multi-dimensional case, we computationally evidence the connection between the agent-based and continuum-based modelling approaches.

@en

We consider a two-dimensional biomorphoelastic model describing post-burn scar contraction. This model describes skin displacement and the development of the effective Eulerian strain in the tissue. Besides these mechanical components, signaling molecules, fibroblasts, myofibroblasts, and collagen also play a significant role in the model. We perform a sensitivity analysis for the independent parameters of the model and focus on the effects on features of the relative surface area and the total strain energy density. We conclude that the most sensitive parameters are the Poisson’s ratio, the equilibrium collagen concentration, the contraction inhibitor constant, and the myofibroblast apoptosis rate. Next to these insights, we perform a sensitivity analysis where the proliferation rates of fibroblasts and myofibroblasts are not the same. The impact of this model adaptation is significant.

@en

Plastic (permanent) deformations were earlier, modeled by a phenomenological model in Peng and Vermolen (Biomech Model Mechanobiol 19(6):2525–2551, 2020). In this manusctipt, we consider a more physics-based formulation that is based on morphoelasticity. We firstly introduce the morphoelasticity approach and investigate the impact of various input variables on the output parameters by sensitivity analysis. A comparison of both model formulations shows that both models give similar computational results. Furthermore, we carry out Monte Carlo simulations of the skin contraction model containing the morphoelasticity approach. Most statistical correlations from the two models are similar, however, the impact of the collagen density on the severeness of contraction is larger for the morphoelasticity model than for the phenomenological model.

@en

The discretisation of the Laplacian results into the well-known Laplace matrix. In the case of a one dimensional problem, an explicit formula for its inverse is derived on the basis of fundamental solutions (Green’s functions) for general boundary conditions. For a linear reaction–diffusion equation, approximations of the inverse are given.@en

Burn injuries can decrease the quality of life of a patient tremendously, because of esthetic reasons and because of contractions that result from them. In severe case, skin contraction takes place at such a large extent that joint mobility of a patient is significantly inhibited. In these cases, one refers to a contracture. In order to predict the evolution of post-wounding skin, several mathematical model frameworks have been set up. These frameworks are based on complicated systems of partial differential equations that need finite element-like discretizations for the approximation of the solution. Since these computational frameworks can be expensive in terms of computation time and resources, we study the applicability of neural networks to reproduce the finite element results. Our neural network is able to simulate the evolution of skin in terms of contraction for over one year. The simulations are based on 25 input parameters that are characteristic for the patient and the injury. One of such input parameters is the stiffness of the skin. The neural network results have yielded an average goodness of fit (R²) of 0.9928 (± 0.0013). Further, a tremendous speed-up of 19354X was obtained with the neural network. We illustrate the applicability by an online medical App that takes into account the age of the patient and the length of the burn.

@en

We consider a mathematical model for skin contraction, which is based on solving a momentum balance under the assumptions of isotropy, homogeneity, Hooke's Law, infinitesimal strain theory and point forces exerted by cells. However, point forces, described by Dirac Delta distributions lead to a singular solution, which in many cases may cause trouble to finite element methods due to a low degree of regularity. Hence, we consider several alternatives to address point forces, that is, whether to treat the region covered by the cells that exert forces as part of the computational domain or as ‘holes’ in the computational domain. The formalisms develop into the immersed boundary approach and the ‘hole’ approach, respectively. Consistency between these approaches is verified in a theoretical setting, but also confirmed computationally. However, the ‘hole’ approach is much more expensive and complicated for its need of mesh adaptation in the case of migrating cells while it increases the numerical accuracy, which makes it hard to adapt to the multi-cell model. Therefore, for multiple cells, we consider the polygon that is used to approximate the boundary of cells that exert contractile forces. It is found that a low degree of polygons, in particular triangular or square shaped cell boundaries, already give acceptable results in engineering precision, so that it is suitable for the situation with a large amount of cells in the computational domain.

@en

A burn wound is a complex systemic disease at multiple levels. Current knowledge of scar formation after burn injury has come from traditional biological and clinical studies. These are normally focused on just a small part of the entire process, which has limited our ability to sufficiently understand the underlying mechanisms and to predict systems behaviour. Scar formation after burn injury is a result of a complex biological system-wound healing. It is a part of a larger whole. In this self-organising system, many components form networks of interactions with each other. These networks of interactions are typically non-linear and change their states dynamically, responding to the environment and showing emergent long-term behaviour. How molecular and cellular data relate to clinical phenomena, especially regarding effective therapies of burn wounds to achieve minimal scarring, is difficult to unravel and comprehend. Complexity science can help bridge this gap by integrating small parts into a larger whole, such that relevant biological mechanisms and data are combined in a computational model to better understand the complexity of the entire biological system. A better understanding of the complex biological system of post-burn scar formation could bring research and treatment regimens to the next level. The aim of this review/position paper is to create more awareness of complexity in scar formation after burn injury by describing the basic principles of complexity science and its potential for burn care professionals. Declaration of interest: The authors gratefully acknowledge financial support from the Dutch Burns Foundation (DBF), Beverwijk, project 19.105. The authors have no conflicts of interest to declare.

@en

We consider a morphoelastic framework that models permanent deformations. The text treats a stability assessment in one dimension and a preservation of symmetry in multiple dimensions. Next, we treat the influence of uncertainty in some of the field variables onto the predicted behaviour of tissue.

@en

We consider a cell-based approach in which the balance of momentum is used to predict the impact of cellular forces on the surrounding tissue. To this extent, the elasticity equation and Dirac Delta distributions are combined. In order to avoid the singularity caused by Dirac Delta distribution, alternative approaches are developed and a Gaussian distribution is used as a smoothed approach. Based on the application that the pulling force is pointing inward the cell, the smoothed particle approach is probed as well. In one dimension, it turns out that the aforementioned three approaches are consistent. For two dimensions, we report a computational consistence between the direct and smoothed particle approach.

@en

The phenomenological model for cell shape deformation and cell migration Chen (BMM 17:1429–1450, 2018), Vermolen and Gefen (BMM 12:301–323, 2012), is extended with the incorporation of cell traction forces and the evolution of cell equilibrium shapes as a result of cell differentiation. Plastic deformations of the extracellular matrix are modelled using morphoelasticity theory. The resulting partial differential differential equations are solved by the use of the finite element method. The paper treats various biological scenarios that entail cell migration and cell shape evolution. The experimental observations in Mak et al. (LC 13:340–348, 2013), where transmigration of cancer cells through narrow apertures is studied, are reproduced using a Monte Carlo framework.

@en

This Special Issue contains a series of fourteen papers that resulted from the Lorentz workshop in Leiden in the period of May 22–25 in 2018.

@en

To deal with permanent deformations and residual stresses, we consider a morphoelastic model for the scar formation as the result of wound healing after a skin trauma. Next to the mechanical components such as strain and displacements, the model accounts for biological constituents such as the concentration of signaling molecules, the cellular densities of fibroblasts and myofibroblasts, and the density of collagen. Here we present stability constraints for the one-dimensional counterpart of this morphoelastic model, for both the continuous and (semi-) discrete problem. We show that the truncation error between these eigenvalues associated with the continuous and semi-discrete problem is of order O(h²). Next we perform numerical validation to these constraints and provide a biological interpretation of the (in)stability. For the mechanical part of the model, the results show the components reach equilibria in a (non) monotonic way, depending on the value of the viscosity. The results show that the parameters of the chemical part of the model need to meet the stability constraint, depending on the decay rate of the signaling molecules, to avoid unrealistic results.

@en