Resilience and repertoire: computational modeling for novel personalized markers

(1)

Bachelor Informatica

Resilience and repertoire:

com-putational modeling for novel

personalized markers

Leonie Rebecca Pelkman

June 15, 2020

Supervisor(s): dr. R. Quax, dr. J. Rector, dr. R. Melis

Inf

orma

tica

—

Universiteit

v

an

Ams

terd

am

(2)

(3)

Abstract

The terms resilience (how well a system is able to cope with perturbations) and repertoire (the set of mechanisms available to cope with those perturbations) are widely used across different fields of science. In gerontology, research into the resilience of a person may aid not only treatment, but also the prevention of declining health in the elderly. Cerebral Autoreg-ulation (CA), responsible for maintaining stable cerebral blood flow (CBF), is a biological system that appears to be maintained in old age, but a decline in CA functioning has been associated with Alzheimer’s disease. Mathematical models for this system have already been developed, but are fairly complicated and highly specific to this one system. Using a first-principles modeling approach, models of increasing complexity were implemented to research whether distributional properties of any system’s response to perturbations hold information about the breadth of its repertoire. It was hypothesized that skewness is inversely associ-ated with the repertoire, since perturbations outside the repertoire are expected to result in larger deviations and longer recovery times of regulated variables. No such association was found in this study, but results indicate that there may be a shift from well-regulated perturbations to perturbations that fall outside the range of repertoire. Around this shift outliers are frequent, causing leptokurtic results which are in some cases positively skewed. Future research should focus on comparing these results to those of more complex models to see if this property holds with greater complexity.

(4)

(5)

Introduction

1.1 Context

1.1.1 Resilience and repertoire

Resilience is a system’s ability to resist and respond to perturbations, response repertoire is defined as a system’s set of mechanisms available to cope with these perturbations. Systems with a broader repertoire at their disposal are expected to have greater resilience to adverse events. Each system has its own set of mechanisms which it uses to counteract fluctuations and return the system to its equilibrium. These definitions are used across a wide range of research fields and thus across a wide range of systems, from very large systems such as in the analysis of the performance of nuclear power plants [6] to processes on the cellular level [5]. One such field is gerontology.

1.1.2 Gerontology

Gerontology is the study of the social, cultural, psychological, cognitive, and biological aspects of aging and is distinguished from geriatrics, which is the branch of medicine that specializes in the treatment of existing disease in older adults.

Aging comes with many changes to physiology. ”Most organ systems show a physiological reduction in function with age, although the rate varies between systems within an individual as well as between individuals. There is reduced redundancy of function and ability to repair” [14]. This narrowing of repertoire causes lower resilience in older people: the older a person gets, the higher the chance becomes that they will die (figure 1.1). In 2017 and 2018 people over 80 years old accounted for less than 5% of the population of the Netherlands [19], but over 55% of the people who died in those years fell in this age group [20].

A decline in physical functioning is another important aspect of aging and causes many people to become reliant on others in old age. In 2019, over 115 thousand people lived in nursing homes in the Netherlands [18]. Healthcare and care for the elderly, and the improvement thereof, require a large budget: In 2017, the Dutch government set aside AC2.1 billion to improve the quality of nursing home care [15] and in 2018, AC340 million was made available for the elderly to keep living independently for longer. Not only does their decreasing independence create a financial burden on society, it can also be a source of stress for them and their caretakers. In that same year, approximately 4.4 million people were caregivers for family, friends or acquaintances in the Netherlands [17]. Although a relatively small proportion of them feels overburdened (10%), this still amounts to 440 thousand people.

Additionally, the population of the Netherlands is aging [17], which means that the challenges we face today will become even greater in the upcoming years. Thus, the narrowing of repertoire and resulting reduction of resilience in the elderly impacts both individuals and society as a whole in multiple significant ways, making it an important subject of study.

(8)

Figure 1.1: Age distribution of deaths in the Netherlands in 2017 and 2018 [20] and the age distribution of residents in the Netherlands in January of 2017 and 2018 [19].

1.1.3 Resilience and repertoire in gerontology

Even within gerontology, resilience and repertoire are characterized and used in many different ways and subfields. Clark et al. (2011) [3] recognizes that in gerontology, there are biological, psychological and social contributions to a person’s resilience and visit several definitions of resilience which have been used in research in the past. These focus mainly on the psychological and social contributions to resilience of a person as a whole. ”However, there has been only limited research in humans on resilience of responses to physiologic or pathological stressors, and clinical and physiologic factors influencing these responses” [9]. This project will therefore be focused on physical resilience.

1.1.4 Physical resilience

Whitson et al. (2016) [24] define physical resilience as ”a characteristic at the whole person level which determines an individual’s ability to resist functional decline or recover physical health following a stressor”. They also propose that the term physiologic reserve be used as the potential capacity of a cell, tissue or organ system to function beyond its basal level in response to alterations in physiologic demands.

As mentioned earlier, aging comes with a progressive decline in physiological functioning of most organ systems, with differences within and between individuals. Navaratnarajah and Jackson (2016) describe age-related changes in physiology in the cardiovascular system, the nervous system, the kidneys, the respiratory system, the gastro-intestinal system, the skin, the haematological system, the endocrine system, the musculoskeletal system, thermoregulation and the immune system [14] [2]. They further state that although these changes don’t necessarily lead to pathology, ”they undoubtedly predispose to pathology” [14]. This statement suits the fact that multimorbidity (or the co-existence within an individual of more than one long-term condition) is common in older adults and is strongly associated with physical functioning, mobility limitations and mortality [23].

For many physiologic systems it is currently impossible to estimate the breadth of the repertoire and quantify the corresponding physiologic reserve. Furthermore, in some if not all -systems, a directional component in the perturbations determines how well a person is able to respond. For example, older people are more vulnerable to falls in a sideways direction than

(9)

falls in a forward or backward direction [13]. Directionality doesn’t need to be literal, it’s simply used to indicate that aspect of the perturbation determine whether it fits within the repertoire. For example: the magnitude of a perturbation may be small enough for the system to handle or too large for mechanisms in the repertoire to return the system to homeostasis. Many factors contribute to physical resilience and its decline is characterized as a heterogeneous process [14]. The directional vulnerability to falls found in the study by Mille et al. (2013) [13] was not present in younger adults. The heterogeneity found in both physical resilience and its decline supports the notion of directionality.

1.1.5 Research needs

Hadley et al. (2017) [9] identified a number of research needs to gain a better understanding of physical resilience, including simulation tests: simulating (within- and between-person level) responses to simulated perturbations can help extract indicators of a narrow response repertoire. Whitson et al. (2011) [24] identify three key areas that warrant future research: measurement of resilience, understanding how system-level reserve contributes to whole-person resilience, and interventions that promote physical resilience in the face of health challenges.

1.1.6 Cerebral Autoregulation

Cerebral autoregulation is a regulatory system that ensures that cerebral blood flow (CBF) remains stable over a wide range of systemic blood pressure (BP) levels [22]. Both static (sCA) and dynamic (dCA) cerebral autoregulation have been a subject of research in the past. ”The classic studies of CA are performed in the semi-steady state” [22], tracking changes in stable blood pressure levels that occur naturally over time and the effect that that has on cerebrovascular resistance. Studies of dynamic CA quantify fast modifications in the velocity of CBF (CBFV) in a major cerebral artery, in relation to rapid alterations in BP within the upper and lower limits of sCA. In terms of resilience and repertoire, alterations in BP within the limits of sCA are perturbations that fall within the repertoire, whereas blood pressure levels that are outside this range do not. CA functions like a high-pass filter in the frequency range of 0.07 to 0.30 Hz; Rapid BP changes are transmitted to CBFV, whereas slow BP changes are filtered [22]. The resilience (stability in the face of perturbations) of CBF depends on its repertoire: cerebral autoregulation.

Although van Beek et al. (2008) concluded that CA is mostly retained in old age, medical conditions that arise from impaired cerebral perfusion (such as syncope and falls caused in part by orthostatic hypotention) are prevalent in the elderly [22]. ”Autoregulation studies can con-tribute valuable clinical information regarding CBF control by identifying or confirming patients at high risk for impaired cerebral perfusion” [22].

Mathematical models that describe dCA have been developed already [21]. In her masters thesis, Jansen [10] implemented a slightly adapted version of this model. Her research was focused on the differences in CA between healthy older adults and older adults with Alzheimer’s disease (AD) and mild cognitive impairment (MCI, a risk factor for AD). One noteworthy finding was that the autoregulatory time constant increases with the severity of the dementia, meaning that it took longer for the CBF to return to normal than in the healthy test subjects. This can be seen as a reduction of resilience and was attributed to the arterial compliance (the arteries’ capability for constriction and dilation) being less capable to respond timely to higher frequency oscillations. The decline in functioning of one of the mechanisms of the repertoire results in a narrower repertoire and thus lower resilience, which brings us to the research question:

1.2 Research question

Non-invasive, physiologic time-series data can be used to simulate physiologic systems, which can be used to make inferences about how these would react to perturbations. For any sys-tem, compared to perturbations within the repertoire, perturbations outside the range of the repertoire are expected to cause larger deviations and require longer recovery times (as seen in

(10)

AD patients in Jansen’s thesis for example [10]), since they require wider and slower secondary and tertiary response mechanisms to counteract the effect. Thus the skewness of the sampled response distribution is hypothesized to be inversely associated with the repertoire. The research question is therefore:

”Do the distributional properties of time-series data in response to perturbations hold information about the breadth of a person’s repertoire, and thus their resilience?” A first-principles approach is used to find an answer for this question. First-principles modeling is an approach that requires no experimental data, but instead attempts to extrapolate and predict patterns using already available knowledge about a system. This is a time- and cost efficient way of determining whether the hypothesis can be confirmed and be used in follow-up research in physical resilience, and has been proven useful in other modeling for other fields, such as the development of new materials [25].

1.3 Ethics

1.3.1 Positive implications

”Physical resilience as a paradigm may offer a next step to take geriatric medicine to a higher level” [8]. Improving the current knowledge of physical resilience may lay the foundations for bet-ter care for the elderly: preemptive screening, improved recovery and new inbet-tervention methods would considerably improve quality of life in older adults.

1.3.2 Important considerations

Although because of time constraints it was not used, non-invasive physiologic time-series data from older adults (collected by dr. Jurgen Claassen in a Radboud University study) was provided to me by my supervisor dr. Jerrald Rector. All data was anonymized before it was entrusted to me, so the privacy of the participants was respected.

Furthermore, one of the goals of researching physical resilience is to improve the ability to measure a person’s resilience. If such measures become widely accepted, they may be used to determine whether they are resilient enough to benefit from certain procedures, such as hip replacement surgery. That means that people who are now excluded from a procedure based on their age, could then undergo them based on their resilience measure. However, it also means that people who would be eligible for a procedure with the current policies, would then be excluded. Since such decisions have a major impact on people, it’s important that the measures are accurate.

Even though this study is too small to directly impact any such policies, it’s important to keep in mind that every study on the road to that goal needs to provide valid and reliable results.

(11)

CHAPTER 2

Theoretical background

2.1 System Dynamics

System dynamics is a computer-aided methodology and mathematical modeling technique, de-signed to approach dynamic systems. Its basic structure consists of a system of coupled differen-tial (or integral) equations. Simulating such systems is accomplished by iterating over discrete intervals of time with length ∆t, and finding the output at each timestep. The input for such a model is typically a value or a vector of values, the output is a value or a vector of values of the same length. System dynamics can be used to simulate behavior of complex dynamic systems.

Dynamical systems are often visualized using block diagrams. Each part of the system can be described by its transfer function which takes an input (an incoming arrow) and generates an output (an outgoing arrow) and is visualized as a block connected to the rest of the system with arrows. Figure 2.1 shows a basic control system where the difference between the target output y0 and the actual output y1 is used as an error, a controller (explained below) then adapts this to a controlled signal, and a plant uses this signal to generate y1. y1 is then subtracted from the output that comes in at the next timestep and the whole process is repeated.

Figure 2.1: Basic control system with a single feedback loop, input y0 and output y1. The general transfer function for a closed-loop system with negative feedback such as this is shown in equation 2.1, with G(s) the transfer function of the controller and the plant (the open-loop parts, and H(s) the transfer function for the feedback loop. In this case the feedback loop doesn’t alter the output of the plant, so H(s) = 1

T (s) = y1(s) y0(s) = G(s) 1 + G(s) · H(s) = G(s) 1 + G(s) (2.1)

2.1.1 Ordinary Differential Equations (ODEs)

An ordinary differential equation is a differential equation containing one or more functions of a single independent variable and the derivatives of those functions, as opposed to a partial

(12)

differential equation, which is derived from one or more functions with multiple independent variables. Solving an ODE means determining what function satisfies, or which functions satisfy the equation. An alternative approach is to use an iterative algorithm which, instead of solving the equation, calculates the slope and adds it to the output at each timestep, as described above. The forward Euler method is such an algorithm. The ODE is provided in the following form:

dy(t)

dt = f (t, y(t)) (2.2)

with an initial value y(t0) = y0. A numeric solution is found by replacing the derivative on the

left-hand side of the equation with a finite difference approximation, then solving for y(t + h): dy(t)

dt ≈

y(t + h) − y(t)

h (2.3)

y(t + h) ≈ y(t) + h ·dy(t)

dt (2.4)

y(n + 1) ≈ y(n) + h ·dy(t)

dt (2.5)

where h is the step size. A smaller step size increases the accuracy, but also the computation cost.

2.1.2 Stochastic Differential Equations (SDEs)

Stochastic differential equations differ from ODEs in one key way: they contain a random com-ponent, making the solution a stochastic process. An ODE can be transformed into an SDE by adding a Wiener process, or Brownian motion.

A Wiener process or Brownian motion is a stochastic variable Wtwhose value changes over time

in an uncertain way. Since the SDEs are solved in discrete time, the Wiener process is also discretized. The change in the value per timestep can be characterized as ∆Wt= Wt+∆t− Wt

and not only contains a random component, but also takes the length of ∆t into account: ∆Wt= εt

√

∆t (2.6)

in which εtis normally distributed random variable with expectation = 0 and variance = 1. This

characterization of ∆Wt implicates that although the changes in the value are independent of

each other, the absolute value of Wt is dependent on the previous timestep. The equation for

WT at the final timestep T is shown in equation 2.7.

WT = (ε0+ ε∆t+ ... + εT −∆t) ·

√

∆t (2.7)

WT is normally distributed with expectance = 0 and variance = T

The Euler-Maruyama method is used to iteratively solve SDEs, it is a generalization of the forward Euler method for ODEs to SDEs, differing from the forward Euler method only in the addition of the stochastic term.

2.1.3 Frequency domain

Many biological systems are complex and to achieve an acceptably reliable model, second-order systems are often used to describe them [4]. Since it can be difficult to solve those increasingly complex transfer functions, the transfer functions are determined for the frequency domain, rather than the time domain. In order to find a function in the frequency domain which corre-sponds to a function in the time domain, a Laplace transform is used. Its inverse is used to find a function in the time domain which corresponds to a function in the frequency domain.

(13)

Laplace transforms are functions whose general form is shown in equation 2.8. Its inverse is shown in equation 2.9. Typically the counterpart of a lowercase function of t is denoted as an upper case function of s, where s is a complex number consisting of both real and imaginary parts. To prevent having to reinvent commonly used Laplace transforms each time they’re needed, Laplace transform tables exist which contain functions in the time domain and their counterparts in the frequency domain. The inverse transforms especially are difficult to find without Laplace transform tables. If necessary, a function F (s) may be algebraically reshaped to resemble one or more functions in the table.

L{f (t)} = F (s) = Z ∞

0

f (t) · e−stdt (2.8)

L−1{F (s)} = f (t) (2.9)

Equation 2.10 shows how the transfer function T (s) for any part of a system can be found using input R(s) and output Y (s).

T (s) = Y (s)

R(s) (2.10)

The general form for a second-order system transfer function in the frequency domain is given in equation 2.11. Here ωndenotes the natural frequency of a process and ζ is known as the damping

factor, which determines the accuracy of the response. ζ < 1 results in an underdampened system in which the system’s response to perturbations is quick, but may also result in larger initial overshoot of the target value. ζ > 1 gives an overdampened system. There is no overshoot but the response is slow. ζ = 1 results in the maximum speed of response for which no overshoot occurs. T (s) = ω 2 n s2_{+ 2ζω} ns + ωsn (2.11)

2.1.4 Control

There are two types of control: loop control and feedback- or closed-loop control. In open-loop control, the control action is determined independently from the output of the process at earlier timesteps. If either the output or some part of the output is used as part of the input (such as in figure 2.1, this is known as feedback. The type of control is then feedback- or closed-loop control. Within closed-loop control there are different types of control mechanisms:

on/off control is a simple form of control which uses a threshold to determine whether an action takes place (on) or not (off).

Proportional control applies a correction to the controlled variable that is proportional to the difference between the measured value and the desired value (the error). Proportional control is characterized by equation 2.12, where u(t) is the output, KP is proportional gain and e(t) is the

error at time t.

u(t) = KP · e(t) (2.12)

To find the transfer function for a proportional controller, Laplace transform is applied to both sides of the equals sign, after which equation 2.10 is applied with input E(s) and output U (s):

L{u(t)} = L{KP· e(t)}

U (s) = KP · E(s) (2.13)

U (s)

(14)

Integral control is proportional to the integral of the error, rather than to the error itself as it is with proportional control. Integral control is characterized by equation 2.14. Again, u(t) is the input, KI is the gain and e(t) is the error at time t.

u(t) = KI·

Z

e(t)dt (2.14)

The transfer function for the integral controller can be found using the same procedure as with the proportional controller: Laplace transform is applied to both sides of the equals sign, after which equation 2.10 is applied with input E(s) and output U (s):

L{u(t)} = L{KI· Z e(t)dt} U (s) = KI· E(s) s (2.15) U (s) E(s) = KI s

Derivative control is proportional to the derivative of the error, rather than to the error itself or its integral as with proportional control and integral control respectively. Derivative control is characterized by equation 2.16. As before, u(t) is the input, KD is the gain and e(t) is the error

at time t.

u(t) = KD·

de(t)

dt (2.16)

The transfer function is found in equation 2.17, using the same procedure as with proportional and integral control.

L{u(t)} = L{KD· de(t) dt } U (s) = KD· s · E(s) (2.17) U (s) E(s)= KD· s

Different types of controllers may be combined to create new kinds of controllers, such as PD, PI or PID controllers. This is done by summing their transfer functions, resulting in a new transfer function for the combined controllers.

2.2 Cerebral Autoregulation

Approximately 20% of the blood flowing from the heart is pumped to the brain. Two sets of arteries provide the brain with blood: the left and right vertebral artery (VA), supplying blood to the caudal portion of the brain (at the back of the head) and the left and right internal carotid artery (ICA), which supplies blood to the rostral portions (on the side of the face) [7] [11]. The ICAs branch off from the common carotid arteries, the bases of which contain stretch receptors that respond to the drop in blood pressure upon standing. The VAs merge into the basilar artery. The left and right ICA and branches of the basilar artery are linked together through the circle of Willis (figure 2.2), which can maintain perfusion of the brain even if blood flow through one part is limited [1]. The circle of Willis is made up of the left and right Anterior Cerebral Arteries, the Anterior Communicating Artery and the Posterior Communicating Arteries. The left and right Middle Cerebral Arteries (MCAs) are also connected to the circle of Willis [16].

The MCAs provide blood for approximately 70% of the brain [10]. Transcranial Doppler ultrasonography (TCD) is a non-invasive measuring tool used to find the flow-velocity in the MCA, denoted as CBFV. It’s assumed that changes in CBFV reflect changes in CBF, making CBFV a useful measure for cerebral autoregulation research [22].

(15)

Figure 2.2: The circle of Willis and connected arteries. This image was taken, without changes, from Anatomy and Physiology by Betts et al. (2013) [1]

As mentioned earlier, Ursino and Lodi [21] created a mathematical model for dynamic cerebral autoregulation, which was implemented by Jansen [10] in her masters thesis. This model was created using a top-down approach, starting by analysis of the system to be modeled, reducing complexity where this is deemed possible without reducing the accuracy of the model too much. Although the model is comparatively simple, it still contains multiple feedback loops and many parameters [10].

Because the model is highly specific to dCA, no knowledge about physiologic systems in general can be gained from it. First-principles modeling, which uses a bottom-up approach, may be more suitable to make general statements about physiologic reserve, and in turn, physical resilience. Comparisons with specific models such as this one can be used to validate first-principles models.

(16)

(17)

CHAPTER 3

Methods

3.1 Generic models of increasing complexity

3.1.1 Self-regulating variable

In order to gain an understanding of system dynamics, I started with a simple ordinary differential equation (ODE) (equation 3.1). Here b is a balance point, KP controls the speed at which y(t)

converges to (with KP > 0) or diverges from (with KP < 0) b. The difference b − y(t) between

the balance point b and the y-value at time t is the error e(t) at time t.

y0(t) = KP· (b − y(t)) = KP · e(t) (3.1)

The numerical solution for this equation is found by applying the forward Euler method, resulting in equation 3.2.

y(t + 1) ≈ y(t) + KP· (b − y(t) · dt = y(t) + KP · e(t) (3.2)

The equation can be changed to a stochastic differential equation (SDE) by adding random noise (equation 3.3). A stochastic Wiener process W is used, its change ∆Wt at time t is added at

each timestep.

y0(t) = KP· (b − y(t)) + ∆Wt= KP · e(t) + ∆Wt (3.3)

Applying the Euler-Maruyama method, the numerical solution for this SDE is given in equation 3.4.

y(t + 1) ≈ y(t) + KP· (b − y(t)) · dt + ∆Wt= y(t) + KP · e(t) · dt + ∆Wt (3.4)

With these equations, the response is proportional to the perturbation regardless of the magni-tude of said perturbation. Every perturbation can therefore be said to be within the repertoire of this system. In section 3.1.3, equations are proposed where the magnitude of a perturba-tion funcperturba-tions as a type of direcperturba-tionality, where the system responds differently based on the magnitude of the perturbation.

3.1.2 Multiple regulating variables

In any given (biological) process, multiple factors can determine its outcome. These factors may have different effects on the process and are therefore represented by multiple separate variables. The ODE in equation 3.5 contains two regulating variables KP 1 and KP 2 which each influence

the remaining error e(t).

(18)

Applying the forward Euler method in the same manner as in section 3.1.1, the numerical solution in equation 3.6 is found.

y(t + 1) ≈ y(t) + KP 1· (b − y(t)) + KP 2· (b − y(t)) (3.6)

≈ y(t) + KP 1· e(t) + KP 2· e(t)

In this example, only one extra regulating variable KP 2 was added, but more variables may be

added in the same manner.

Adding random noise in the form of a Wiener process, and finding the numerical solution for the resulting SDE is also done in the same way as in section 3.1.1, and shown in equations 3.7 and 3.8.

y0(t) = KP 1· (b − y(t)) + KP 2· (b − y(t)) + ∆W = KP 1· e(t) + KP 2· e(t) + ∆W (3.7)

y(t + 1) ≈ y(t) + KP 1· (b − y(t)) + KP 2· (b − y(t)) + ∆W (3.8)

≈ y(t) + KP 1· e(t) + KP 2· e(t) + ∆W

In these equations, both regulating variables remain proportional to the perturbation regardless of the magnitude of said perturbation. They could therefore be simplified by summing KP 1 and

KP 2 and using the result as KP in equations 3.1 and 3.3. So as with the single self-regulating

variable, each perturbation - regardless of its magnitude - falls within the repertoire of the system. In section 3.1.3 equations are proposed where KP 1 and KP 2 have different magnitudes

of perturbations for which they respond optimally.

3.1.3 Directionality

In 1959, Lassen [12] was the first to show an autoregulation curve. A plateau of unchanging CBF regardless of the systemic blood pressure (BP). Outside of the bounds of this plateau however, CBF was shown to be highly dependent on BP. This indicates that there is a range of BP values for which CA functions adequately and that the magnitude of the autoregulatory response has a maximum. Implementing such a response (magnitude as directionality) can be done in several ways and will be discussed separately for a single self-regulating variable and for multiple regulating variables.

Directionality with a single self-regulating variable. Three methods of using a maximum response magnitude were implemented:

• Using the previously implemented equation 3.4 but manually choosing a maximum value for KP· e(t) · dt.

• Using equation 3.4, but with a normally distributed KP depending on the error size.

• Using a sigmoidal function as the response.

The implementations are described in more detail in the setup of the experiments.

Directionality with multiple regulating variables was implemented using a maximum response magnitude in the same way as for the single self-regulating variable. The implementations are described in more detail in the setup of the experiments.

(19)

CHAPTER 4

Experiments

4.1 Setup

All models were written in Python using Jupyter notebooks, which allows for the models and their results to be shown in a structured manner.

4.1.1 Self-regulating variable

For the first three experiments, equation 3.4 is implemented. This function takes six arguments: • y0: The value of y at t = 0.

• KP

• b

• n: The amount of timesteps. • dt: The length of each timestep.

• sims: This variable determines how often the equation is executed.

Since the modeled equation is stochastic, the output differs each time the equation is executed. Since the change in a Wiener process ∆W is dependent on the size of dt, the same value for dt is used in each experiment.

The function returns an array of size (sims,) containing y(T ) for each time the simulation is executed. The experiments described below are run with n = 1000, dt = 0.02 and sims = 10000. The chosen values for n and dt result in a final time T = n · dt = 20.

The distribution of the Wiener processes at t = T is visualized first. This is achieved by setting KP to 0 and using y0 = 0 and b = 0. Because of its definition in equation 2.7, the results are

expected to be normally distributed with a mean of 0 and a variance of 20.

The influence of KP on the distribution of y at t = T is shown next. Again using y0 = 0 and

b = 0, the same experiment is performed with different values for KP.

In section 1.2 the skewness of the sampled response distribution is hypothesized to be inversely associated with the repertoire, and since each perturbation - regardless of its magnitude - falls within the repertoire as explained in section 3.1.1, the distributions are hypothesized to be normally distributed and thus have negligible skewness. The variance is expected to become lower as KP increases, since a larger KP means a stronger reduction of the deviation from

balance b at each timestep.

However, should KP become so large that KP · dt > 1, the error e(t) will be multiplied by

(20)

KP · dt > 1, variance is expected to increase as KP is increased. If KP· dt = 2 the overshoot is

expected to be equal to e(t) before ∆Wtis added, resulting in a new error e(t+1) = −e(t)+∆Wt.

In this case the absolute value of the error (and thereby also the variance of y(T )) is expected to be no different than with KP = 0. Finally, if KP · dt > 2, the overshoot is larger than the

original error and a larger variance compared to KP = 0 is expected.

The influence of a large perturbation on the distribution of y at t = T is shown by repeating the experiment using KP = 0.1 and b = 0 with different values for y0. Establishing what will be

characterized as a large perturbation is needed in order to determine which values to use for y0. The previous experiments all contained the change in the Wiener process ∆Wtwith its variance

of dt = 0.02 at each timestep as perturbations, so a large perturbation should be significantly larger than that.

As in the previous experiment, the magnitude of the perturbation is hypothesized to be irrelevant for the skewness of the response distribution in this equation, since each perturbation falls within the repertoire. The variance is expected to remain approximately equal for each y0, since KP is unchanged. Although all perturbations fall within the repertoire, very large ones may

still require more time for the system to return to normal y-values than given in the experiment. Thus the mean is expected to increase as the tested perturbations become larger.

4.1.2 Multiple regulating variables

For the next experiment, equation 3.8 is implemented. This function takes seven arguments, differing only from the arguments in the implementation of equation 3.4 in that it takes two regulating variables KP 1 and KP 2 instead of just one, KP. The same values for n, dt and sims

are used.

For the next experiments involving directionality, it will be assumed that with KP equalling

KP 1+ KP 2, equation 3.8 yields the same results as 3.4. The goal of this experiment is to show

that this is indeed the case. This is done by choosing a number of values for KP 1 and KP 2

and comparing the resulting distributions of both equations. It’s expected that all results are normally distributed with a variance approximately equal to the variance found in the single self-regulating variable implementation.

4.1.3 Directionality

Directionality with a single self-regulating variable Implementation 1, maximum response magnitude:

In the first implementation, the absolute value of KP · e(t) · dt is compared to a manually set

maximum response magnitude. The smallest is chosen as the absolute value of the response. Whether the response is positive or negative (for non-zero error values) is stored by dividing the error by its absolute value (resulting in either 1 or −1) and multiplying this by the resulting response magnitude. The manual maximum response magnitude is determined while keeping in mind that the Wiener process produces normal perturbations which should always fall within the repertoire and should therefore be regulated normally. A slightly modified implementation of equation 3.4 will be used to determine the largest absolute value for ∆Wtand choose a maximum

response magnitude larger than that.

It is expected that large perturbations now fall outside the repertoire and will result in skewed response distributions, in accordance with the hypothesis that the skewness of the response is inversely associated with the repertoire.

Implementation 2, normally distributed KP:

The second implementation uses a normally distributed value of KP, which is dependent on e(t).

Instead of a value for KP, this implementation takes a mean, standard deviation and surface

for the distribution of KP as input arguments. The mean value corresponds to the error value

for which KP will respond optimally. A second function Kp val is called at each timestep and

calculates KP using the normal probability density function and multiplying the outcome by the

(21)

The expectation is that values for e(t) close to the mean of KP will be well-regulated, but

that once the difference between them becomes larger, KP will be so low that the distribution

will approximate the Wiener process distribution. Implementation 3, sigmoidal function for the response:

In the previous implementation, when the difference between e(t) and the mean of KP becomes

too large, KP’s influence on e(t) is expected to become so small that it’s negligible. Since it

would be counter-intuitive for a system’s response to become smaller as a perturbation increases, another way of determining the response magnitude is considered. A sigmoidal function is used as the response magnitude, rather than using it as KP and multiplying it by e(t) · dt.

A cumulative distribution function of a normal distribution, when plotted, has an s-shape, hence the name sigmoid. Starting from zero, the magnitude of the response increases slowly at first, then quicker as the mean is approached, and finally slowing down as it gets closer to the maximum value. This maximum value (determined by the surface) determines which perturbation magnitudes are classified as being within the repertoire and which ones are classified as not being in the repertoire.

It’s expected that this implementation will yield results similar to those from implementation 1, where the response magnitude is also bounded by a maximum, but does not decrease after a certain distance from the balance point is reached like it does in implementation 2.

Directionality with multiple regulating variables Implementation 1, maximum response magnitudes:

Rather than having just one maximum response magnitude, different maxima are used for the separate regulating variables KP 1 and KP 2. After determining the responses and its direction,

equation 3.8 is followed.

Depending on the values for the maxima, the breadth of the repertoire can be adjusted. It is expected that the main difference between the distributions of the implementation with just one regulating variable, and this implementation with multiple regulating variables lies in the error values larger than the maximum for the first regulating variable, but are still smaller than one or more maxima of other regulating variables.

Implementation 2, normally distributed KP 1 and KP 2:

KP 1 and KP 2 are normally distributed, each with their own mean, standard deviation and

surface. After determining the responses and their direction, equation 3.8 is followed.

The expectation is that if there is a valley between the peaks of the normal distributions, this translates to the distribution of the error values and may result in two peaks in the distribution of y(T ).

Implementation 3, sigmoidal functions for KP 1 and KP 2:

KP 1and KP 2are determined using the cumulative distribution function of a normal distribution,

after which equation 3.8 is followed.

It is expected that if there is a valley between the distributions, this may result in two separate peaks in the distribution of y(T ). Another expectation is that the behavior resulting from this implementation will be similar to that from implementation 1, because after reaching its maximum response value, it doesn’t decrease as it does in implementation 2.

4.2 Results

To prevent ambiguity, it should be mentioned that in this project, Fisher’s definition of kurtosis (Pearson’s definition minus 3) is used for kurtosis. Furthermore, a distribution is considered to be normally distributed unless its p-value is smaller than α, with α = 10−3

(22)

4.2.1 Self-regulating variable

The distribution of the Wiener processes at t = T is shown in figure 4.1a. As expected, the Wiener process is normally distributed (z ≈ 0.897, p ≈ 0.639) with mean ≈ 0.025 and variance ≈ 19.466.

The influence of KP on the distribution of y at t = T can be seen in figure 4.1b with KP = 0.1.

From the image it’s already visible that the variance is much lower with KP = 0.1 than with

just the Wiener process, where KP = 0.

Ten values of KP were tested, among which the values for which KP· dt = 1 and KP· dt = 2.

With dt = 0.02, these are KP = 50 and KP = 100 respectively.

Five values of KP were tested for which 0 < KP · dt ≤ 1, namely 0.1, 1, 10, 30 and 50.

The distribution of y(T ) for KP = 0.1 is shown in the histogram (figure 4.1b), from which it’s

already clear that the variance is lower than with just the Wiener processes, for which KP = 0.

As expected all distributions were normal and with KP increasing, the variance decreased from

approximately 5.017 for KP = 0.1, to 0.505, to 0.056, to 0.024 and finally to 0.019 for KP = 50.

Four values of KP were tested for which 1 < KP · dt ≤ 2, namely 60, 70, 80 and 100. Figure

4.1c shows the distribution of y(T ) for KP = 100, for which KP · dt = 2. As expected all

distributions were normal and with KP increasing, the variance increased from approximately

0.020 for KP = 60, to 0.023, to 0.032 and finally to 20.365 for KP = 100, which approximates

the expected variance of unregulated Wiener processes.

Lastly, one value of KP was tested for which KP· dt > 2, namely 101. This also resulted in a

normal distribution, but the variance was extremely high (≈ 7.5 · 1016_{) since for each timestep,}

the new error e(t + 1) equals −101e(t) · dt + ∆Wt= −e(t) + −e(t) · dt + ∆Wt.

The influence of a large perturbation on the distribution of y at t = T is shown in figure 4.1d with y0 = 50.

Five values of y0 were chosen for testing, which all had to be much larger than the per-turbations caused by the Wiener processes. y0 = 10, 50, 100, 200 and 500 were used in this experiment. All distributions were normal, with the variances staying approximately the same. As expected, the mean increased from approximately 1.372 with y0 = 10, to 6.760, to 13.509, to 27.072 and finally to 67.627 with y0 = 500.

4.2.2 Multiple regulating variables

The first values I chose for KP 1 and KP 2 were 0 and 0 (with y0 = 0 and b = 0) to show the

distribution of unregulated Wiener processes as in the self-regulating variable. As expected the result was normally distributed (z ≈ 0.004, p ≈ 0.998) with mean ≈ 0.012 and variance ≈ 19.699. Next, KP 2 = 0 was maintained and several different values for KP 1were used, again with y0 = 0

and b = 0. I chose KP 1= 0.1, 10 and 80 since these were already tested in the implementation of

the single self-regulating variable, with two values where KP 1· dt < 1 and one where KP 1· dt > 1.

The same experiment was performed maintaining KP 1= 0, with KP 2= 0.1, 10 and 80. All

results were normally distributed. Table 4.1 shows the resulting variances, the slight differences within each value for the regulating variable can be explained by the stochastic nature of the implemented equations.

Then, nonzero values were selected for KP 1 and KP 2, for which the sums equal the values

for KP which were used before. For KP 1+ KP 2, 0.1 = 0.02 + 0.08, 10 = 8 + 2 and 80 = 20 + 60

were used. All results were normally distributed, their variances are indeed similar to the results from all previous experiments and are shown in the bottom row of table 4.1.

4.2.3 Directionality

As in the previous experiments, b = 0, n = 1000, dt = 0.02 and sims = 10000 were used. The same values of y0 were used for testing in each following implementation: y0 = 0, 10, 50, 100, 200 and 500.

(23)

(a) The distribution of 10000 Wiener processes at time T with n = 1000 and dt = 0.02.

(b) The distribution of 10000 values for y(T ) with n = 1000, dt = 0.02, KP = 0.1 and y0 = 0.

(c) The distribution of 10000 values for y(T ) with n = 1000, dt = 0.02, KP = 100 and y0 = 0.

(d) The distribution of 10000 values for y(T ) with n = 1000, dt = 0.02, KP = 0.1 and y0 = 50

Figure 4.1: The distribution of unregulated Wiener processes shows higher variance than the distribution of y(T ) with KP = 0.1. For KP · dt = 2, which is the case in subfigure c, the

variance is approximately equal to unregulated Wiener processes. The distribution for KP = 0.1

maintains approximately equal variance with different values for y0. Regulating variable 0.1 10 80

KP 5.017 0.056 0.032

KP 1, (KP 2= 0) 4.854 0.055 0.032

KP 2, (KP 1= 0) 4.926 0.057 0.031

KP 1+ KP 2 4.995 0.055 0.031

Table 4.1: The variances of the distributions of y(T ) for several values of KP are similar to those

resulting from KP 1+ KP 2.

Directionality with a single self-regulating variable Implementation 1, maximum response magnitude:

The maximum absolute value for ∆Wtwas 0.726, so the maximum response magnitude was set

at 1. KP was set at 50 so that KP· dt = 1.

Figure 4.2a shows the distribution for y0 = 0. Contrary to the expectation, the tests with the largest perturbations were all normally distributed. The tests for y0 = 0 and y0 = 10 were not, the distribution of y0 = 0 being negatively skewed (z ≈ −5.490, p ≈ 4.023 · 10−8) and leptokurtic (z ≈ 18.311, p ≈ 6.724 · 10−75) and the distribution of y0 = 10 skewed in the other direction (z ≈ 9.568, p ≈ 1.089) and also leptokurtic (z ≈ 21.960, p ≈ 6.950 · 10−107).

Implementation 2, normally distributed KP:

The mean for KP was chosen to be 0, since it’s expected that regulatory mechanisms are most

effective around the balance value. For the standard deviation, the value 1 was selected because all absolute values for ∆Wtare expected to be lower than that. A surface value of 125 was used,

so that when e(t) = 0, KP ≈ 50 and thus KP · dt ≈ 1. Because the mean equals 0 and the

distribution is symmetrical, it’s not necessary to use the absolute error, so no measures need to be taken to remember whether the error is positive or negative.

The distribution with y0 = 0 had a variance of approximately 0.020, comparable to the vari-ance with KP = 50 in section 4.2.1. Figure 4.2b shows the results for y0 = 10. The two peaks

(24)

visual inspection seems to confirm the hypothesis that the values close to the mean of KP are

regulated, whereas for other values the distribution is more similar to that of unregulated Wiener processes. This is also corroborated by the properties of the distributions for the remaining val-ues of y0, which are all normally distributed with a variance of approximately 20, their means approximately the same as their corresponding y0-value.

(a) The distribution of 10000 values for y(T ) with n = 1000, y0 = 0 and dt = 0.02, using the first approach to directionality with a maximum response magnitude.

(b) The distribution of 10000 values

for y(T ) with n = 1000, y0 = 10 and dt = 0.02, using the second approach to directionality with a normally distributed KP. µKP = 0,

σKP = 1 and surface = 125.

Figure 4.2: Resulting distributions from implementation 1 and implementation 2 of directionality with a single self-regulating variable.

Implementation 3: sigmoidal function for the response:

As in implementation 1, the maximum response magnitude was set at 1. The mean (0.5) and standard deviation (0.4) were selected by plotting how much over- or undershoot would occur for different error values, the plot for the chosen mean and standard deviation is shown in figure 4.3.

The distribution with y0 = 0 was leptokurtic (z ≈ 16.438, p ≈ 1.018 · 10−60), the distribution for y = 10 was also leptokurtic (z ≈ 28.179, p ≈ 1.062·10−174) and positively skewed (z ≈ 15.565, p ≈ 1.268 · 10−54_{). The remaining distributions were all normally distributed with a variance of}

approximately 20.

These results show some similarities with those of implementation 1, which was expected, but since the results from implementation 1 differed from the expectations, so do these results.

Figure 4.3: The undershoot for small error values with one sigmoidal regulating variable KP,

which has mean = 0.5, standard deviation = 0.4 and surface = 1. Negative values indicate overshoot.

(25)

Directionality with multiple regulating variables Implementation 1, maximum response magnitudes:

The maximum response magnitude for KP 1· e(t) · dt was maintained at 1. KP 2· e(t) · dt =

2 was chosen as a second maximum. This results in a broader repertoire than that used in implementation 1, which used KP · e(t) · dt = 1. KP 1 and KP 2 were both set at 25 so that

(KP 1+ KP 2) · e(t) · dt = 1 for error values smaller than KP 1’s maximum response magnitude.

The distribution for y0 = 0 was normal, but the p-value (0.005) was very close to α, showing a trend towards leptokurticity. y0 = 10 resulted in a normal distribution, but y0 = 50 was positively skewed (z ≈ 20.320, p ≈ 8.621 · 10−92) and leptokurtic (z ≈ 29.499, p ≈ 2.977 · 10−191). The remaining distributions were all normal with a variance approximating 20, corresponding to the distribution of unregulated Wiener processes, but with lower means than the starting value, so shifted.

Implementation 2, normally distributed KP 1 and KP 2:

The normal distribution of KP 1 is characterized by mean = 0, standard deviation = 1 and

a surface of 125, the same distribution as used for the experiment with a single self-regulating variable KP. The distribution of KP 2has a mean of 5, standard deviation = 2 and surface = 125

and is therefore flatter than KP 1. There is a valley between the two peaks of the regulating

variables, where regulation is less effective than at either peak.

Contrary to the expectations, the results were normally distributed for each tested value of y0. The variance was very low for y0 = 0 (variance ≈ 0.020) and y0 = 10 (variance ≈ 0.020) and approximated 20 for all other values of y0, their means approximately the same as their starting value y0.

Implementation 3, sigmoidal functions for KP 1 and KP 2:

The mean and standard deviation of KP 1 were kept from the implementation with a single

self-regulating variable (mean = 0.5, standard deviation = 0.4). The repertoire was broadened by adding KP 2 with mean = 5, standard deviation = 1.5 and surface of 6.5, essentially setting the

maximum response magnitude to 7.5. There is a gap between the lowest undershoot caused by KP 1 and the lowest undershoot caused by KP 2. Figure 4.4 shows the undershoot with these

distributions of KP 1 and KP 2 for error sizes up to 10.

Except the distributions from y0 = 200 and y0 = 500, all distributions were leptokurtic with variances of approximately 0.64. The distributions of y0 = 200 and y0 = 500 were normally distributed with a variance of approximately 20, like unregulated Wiener processes, but with lower means than the starting values, so shifted.

Figure 4.4: The undershoot for error values up to 10 with two sigmoidal regulating variables, KP 1

(mean = 0.5, standard deviation = 0.4, surface = 1) and KP 2 (mean = 5, standard deviation

(26)

4.2.4 Increasing the sample size

In the previous section, the result for implementation 1 of directionality with a self-regulating variable with y0 = 0 yielded a statistically significant skewed distribution, whereas a normal distribution was expected since there was no initial error. The result for y0 = 10 was skewed in the other direction. This indicates that the sample size (sims) may have been too small. Therefore, all experiments from section 4.2.3 were repeated with sims = 20000, doubling the sample size.

Most of the outcomes were replicated and indeed, the negative skewness for y0 = 0 in implemen-tation 1 of directionality with a self-regulating variable was no longer present with the increased sample size. Furthermore, in implementation 1 of directionality with multiple regulating vari-ables, the trend towards leptokurticity with y0 = 0 disappeared.

(27)

CHAPTER 5

Discussion

The aim of this research was to gain insight in the mechanics of autoregulatory systems and to see whether distributional properties hold information on the repertoire of such systems. Below, the results of the experiments and their implications will be examined and limitations of this study and suggestions for future research will be discussed.

5.1 Results

The sampled responses for perturbations within the range of the repertoire were hypothesized to be normally distributed. In the first experiments, before directionality was added and the proportion of the error that was regulated was therefore not dependent on the magnitude of the perturbation, this was indeed the case.

The results from the experiments with directionality are more difficult to interpret. Few of the results were skewed, contrary to the expectations. The results that were skewed occurred with values of y0 between those resulting in distributions with small variance and those that were normally distributed with variance equal to that of Wiener process distributions.

No predictions for kurtosis were made initially, but leptokurticity was measured relatively often in implementations 1 and 3 for both a single self-regulating variable and for multiple regulating variables. When leptokurticity occurred for the single self-regulating variable and in implementation 3 for multiple regulating variables, it was always in the lowest values of y0 that were tested. However, in implementation 1 for multiple regulating variables, leptokurticity was only observed with y0 = 50.

One possible explanation for these results is that both skewness and leptokurticity occur around a shift from well-regulated perturbations within the repertoire, to less-regulated pertur-bations outside the range of the response repertoire. This would mean that a large amount of outliers can be observed in a wide range around this shift and that in a smaller range around it, these outliers occur more in the direction of the error than the other way, thus explaining the positively skewed distributions.

The normally distributed results for high values of y0 across all three implementations with a maximum response magnitude have a straightforward explanation. The error values for these perturbations all fall outside the repertoire (i.e. they are larger than the maximum response magnitude), even after T timesteps. Therefore, at each timestep and in every simulation the same value is subtracted from the error, essentially changing nothing about the distribution of the Wiener process. In the case of implementation 2, where falling outside the repertoire entails that no correction is applied at all, the mean remains approximately equal to y0. In implementations 1 and 3, the uniform correction at each timestep results in the mean being shifted closer to the balance point.

(28)

5.2 Limitations

A pitfall of this study is that it’s difficult to validate the results obtained in this study. Since the same experiments haven’t been performed on experimental data, no comparison could be made to see if any such distributional properties occur in more complex models as well.

Also, I hypothesized the normal distribution of perturbations within the repertoire. Although this was confirmed, it’s unclear if the cause for this distribution was indeed the fact that the perturbations were within the repertoire, or if some other cause hasn’t been considered, since many distributions which fell outside the repertoire were also normally distributed.

5.3 Conclusion and future work

Summarizing, although the distributions were not skewed for perturbations outside the reper-toire as hypothesized, a different effect was found. The leptokurtic and - in some cases - skewed results may indicate that there may be a shift between well-regulated and unregulated pertur-bations around which more outliers occur, which are either on both sides, or occurring more in the direction of the perturbation.

Future studies should look into this shift and whether this property holds for models with greater complexity.

Secondly, it would be interesting to see what the distributions look like for the amount of time it takes for different perturbations to return to a value below an error threshold.

Finally, I received the python code of Jansen’s master thesis [10], courtesy of Swetta Jansen herself, and physiologic time series data of older adults obtained by dr. Jurgen Claassen in a Radboud study, provided to me by my supervisor dr. Jerrald Rector. Unfortunately, because of time constraints I was unable to adapt the code for an added perturbation, to perform ex-periments with this code and compare the results to the results from my own models. This comparison could play a big part in validating the results from this study.

(29)

Bibliography

[1] J. G. Betts et al. Anatomy and Physiology. Houston, Texas, USA: OpenStax, Apr. 2013. url: https://openstax.org/books/anatomy-and-physiology/pages/1-introduction. [2] Marcia A. Blackman and David L. Woodland. “The narrowing of the CD8 T cell repertoire

in old age”. In: Current opinion in immunology 23.4 (2011), pp. 537–542.

[3] Phillip G. Clark et al. “What do we know about resilience in older adults? An exploration of some facts, factors, and facets”. In: Resilience in Aging. Springer, 2011, pp. 51–66. [4] Carlo Cosentino and Declan Bates. An Introduction to Feedback control in systems biology.

Crc Press, 2011.

[5] Mathias Faure et al. “Tolerance to maternal immunoglobulins: resilience of the specific T cell repertoire in spite of long-lasting perturbations”. In: The Journal of Immunology 163.12 (1999), pp. 6511–6519.

[6] Dominic Furniss et al. “A resilience markers framework for small teams”. In: Reliability Engineering & System Safety 96.1 (2011), pp. 2–10.

[7] Michael S. Gazzaniga, Richard B. Ivry, and George R. Mangun. Cognitive Neuroscience, The Biology of the Mind. 3rd ed. WW Norton & Company, 2009.

[8] Sanne M. W. Gijzel et al. “Resilience in clinical care: Getting a grip on the recovery potential of older adults”. In: Journal of the American Geriatrics Society 67.12 (2019), pp. 2650–2657.

[9] Evan C. Hadley, George A. Kuchel, and Anne B. Newman. “Report: NIA workshop on measures of physiologic resiliencies in human aging”. In: The Journals of Gerontology: Series A 72.7 (2017), pp. 980–990.

[10] Swetta Anne Jansen. “Parameter Optimization of a Mathematical Model Simulating Cere-bral Autoregulation Using Data of Patients with MCI and AD”. Universiteit van Amster-dam Vrije Universiteit AmsterAmster-dam, June 2019.

[11] Payne Stephen John. Cerebral Blood Flow and Metabolism: A Quantitative Approach. World Scientific, 2017.

[12] Niels A Lassen. “Cerebral blood flow and oxygen consumption in man”. In: Physiological reviews 39.2 (1959), pp. 183–238.

[13] Marie-Laure Mille et al. “One step, two steps, three steps more. . . directional vulnera-bility to falls in community-dwelling older people”. In: Journals of Gerontology Series A: Biomedical Sciences and Medical Sciences 68.12 (2013), pp. 1540–1548.

[14] Arunraj Navaratnarajah and Stephen H. D. Jackson. “The physiology of ageing”. In: Medicine 45.1 (2016), pp. 6–10.

[15] _{Rijksoverheid Nederland. Zorg voor ouderen in verpleeghuizen verbeteren. 2017. url: https:} / / www . rijksoverheid . nl / onderwerpen / verpleeghuizen - en - zorginstellingen / zorg-ouderen-verpleeghuizen-verbeteren.

[16] Stephen Payne. Cerebral Autoregulation: Control of Blood Flow in the Brain. Springer, 2016.

(30)

[17] Programma Langer Thuis. Tech. rep. Ministerie van Volksgezondheid, Welzijn en Sport, 2018.

[18] Centraal Bureau voor de Statistiek (CBS). 115 duizend mensen in verzorgings- of ver-pleeghuis. 2020. url: https://www.cbs.nl/nl-nl/achtergrond/2020/13/115-duizend-mensen-in-verzorgings-of-verpleeghuis.

[19] Centraal Bureau voor de Statistiek (CBS). Bevolking; geslacht, leeftijd en burgerlijke staat, 1 januari. 2019. url: https : / / opendata . cbs . nl / statline / ?dl = 308BE # /CBS / nl / dataset/7461bev/table.

[20] Centraal Bureau voor de Statistiek (CBS). Bevolking; geslacht, leeftijd en burgerlijke staat, 1 januari. 2019. url: https : / / opendata . cbs . nl / statline / # / CBS / nl / dataset / 70895ned/table?fromstatweb.

[21] Mauro Ursino and Carlo Alberto Lodi. “A simple mathematical model of the interaction between intracranial pressure and cerebral hemodynamics”. In: Journal of Applied Physi-ology 82.4 (1997), pp. 1256–1269.

[22] Arenda H. E. A. Van Beek et al. “Cerebral autoregulation: an overview of current concepts and methodology with special focus on the elderly”. In: Journal of Cerebral Blood Flow & Metabolism 28.6 (2008), pp. 1071–1085.

[23] Melissa Y. Wei et al. “Physical functioning decline and mortality in older adults with multimorbidity: joint modeling of longitudinal and survival data”. In: The Journals of Gerontology: Series A 74.2 (2019), pp. 226–232.

[24] Heather E. Whitson et al. “Physical resilience in older adults: systematic review and de-velopment of an emerging construct”. In: Journals of Gerontology Series A: Biomedical Sciences and Medical Sciences 71.4 (2016), pp. 489–495.

[25] Yue Yuan, Huabei You, and Luis Ricardez-Sandoval. “Recent advances on first-principles modeling for the design of materials in CO2 capture technologies”. In: Chinese Journal of Chemical Engineering 27.7 (2019), pp. 1554–1565.

Resilience and repertoire: computational modeling for novel personalized markers

Bachelor Informatica