**3. Fitting the simulation model to observational data**

As described in the "Introduction" section, cancer genome analysis demonstrated intratumor heterogeneity and branching evolution of cancer; paticulally, an approach known as multiregion sequencing has been popularly employed for analyzing solid tumors. Here, we introduce a concrete example of a multiregion

sequencing study and explain the utility of cancer evolution simulation when combined with multiregion sequencing data.

In multiregion sequencing, multiple samples obtained from physically separate regions within the tumor of a single patient are analyzed (**Figure 3A**), with two categories of somatic single-nucleotide mutations identified: "founder" and "progressor" mutations (**Figure 3B**). Founder mutations are defined as present in all regions, whereas progressor mutations are defined as present in some regions (note that they are also referred to using different terms in different studies, e.g., public/private or trunk/branch mutations). Founder mutations are thought to accumulate during the early phases of cancer evolution. The common ancestor clone acquires all founder mutations, and then branches into subclones, which accumulate progressor mutations and contribute to forming ITH. Through these multiregion mutational profiles, we can infer an evolutionary history of the cancer by constructing a phylogenetic tree (**Figure 3C**).

As a pioneering study, Gerlinger et al. [17] performed multiregion sequencing, revealing extensive ITH and clonal branching evolution in renal cancer. Ther also identified not only founder mutations in some known driver genes such as *VHL*, but also progressor mutations in other known driver genes such as *SETD2* and *BAP1*. Interestingly, in some cases, different mutations in the same driver gene or genes with the same function were acquired independently. This phenomenon known as parallel evolution also indicates that part of the ITH was generated by Darwinian selection.

Uchi et al. [5] also investigated ITH in nine cases of surgically resected late-stage colorectal tumors by multiregion sequencing to identify founder and progressor mutations in each case. **Figure 4** shows the results obtained from one of the nine cases, which contains 20 samples from the primary lesion and one sample from the metastatic lesion. Note that the progressor mutations showed a mutational pattern that was geographically correlated with the sampling locations. Moreover, they found that mutation allele frequencies, which can be approximately regarded as the proportion of cells with mutations in each region, tended to be lower for progressor mutations than for founder mutations. This observation suggests that the founder mutations existed in all the cancer cells while the progressor mutations existed in only a fraction of the cancer cells in each region. Thus, even in each region, extensive ITH may have existed, which was not captured by the resolution of multiregion sequencing. In addition, most mutations in known driver genes such as *APC* and *KRAS* were identified as founder mutations. However, progressor mutations contain few driver mutations and parallel evolution was not confirmed, which contrasts to the findings obtained in renal cancer. These observations suggest that apart from

#### **Figure 3.**

*Multiregion sequencing (A) DNA samples from multiple regions of a single tumor are analyzed by nextgeneration sequencing. (B) Through multiregion mutation profiling, founder and progressor mutations are identified as common mutations in all regions tested and only restricted regions, respectively. (C) In a phylogenetic tree constructed from the multiregion mutation profile, the trunk and branches correspond to the founder and progressor mutations, respectively. This image originally appeared in [6].*

*Agent-Based Modeling and Analysis of Cancer Evolution DOI: http://dx.doi.org/10.5772/intechopen.100140*

#### **Figure 4.**

*Multiregion sequencing of colorectal cancer. (A) Schema of the tumor subjected to multiregion sequencing. (B) Multiregion mutation profile. The depth of red represents the mutant allele frequency, whereas the colors of the sample labels were prepared so that the similarities of colors represent those of mutation patterns. (C) Phylogenetic tree constructed from the multiregion mutation profile. The time when mutations in known driver genes of colorectal cancer were acquired is indicated along the tree. This image was obtained by modifying a figure that originally appeared in [5].*

Darwinian selection, there are other evolutionary principles generating ITH. To identify these principles, they developed the BEP model as described above; by fitting the BEP model to the multiregion sequencing data, they evaluated the evolutionary principles generating ITH in colorectal cancer.

To fit the simulation model to the observational data, we can employ ABC [18], which constitutes a class of computational methods rooted in Bayesian statistics that can be used to estimate the posterior distributions of model parameters. A common incarnation of Bayes' theorem relates the conditional probability of a specific parameter value *θ* given data *D* to the probability of *D* given *θ* by the rule, *p*ð Þ *θ*j*D* ∝*p D*ð Þ j*θ p*ð Þ*θ* , where *p*ð Þ *θ*j*D* denotes the posterior, *p D*ð Þ j*θ* the likelihood, and *p*ð Þ*θ* the prior. The prior represents beliefs or knowledge about *θ* before *D* is available. To obtain the the posterior, the likelihood function is required. For simple models, an analytical formula for the likelihood function can typically be derived. However, for more complex models, an analytical formula may be elusive or the likelihood function may be computationally very costly to evaluate. Agent-based models also fall into the latter case. ABC methods bypass evaluation of the likelihood function by using summary statistics and simulations, which widen the realm of models for which statistical inference can be considered. ABC has rapidly gained popularity over the last few years, for analizing complex problems arising in biological sciences, e.g., in population genetics, ecology, epidemiology, and systems biology.

In the basic form of ABC, which is known as rejection sampling, we first sample a parameter value (or a combination of parameter values, if there is more than one parameter) from a prescribed prior distribution of the parameter value. Simulated data are then generated from the sampled parameter value. The similarity between the simulated and observational data is evaluated using summary statistics (typically multiple), which is designed to represent the maximum amount of information in the simplest possible form. If the distance of the summary statistics between the simulated and observational data is below a tolerance parameter, the parameter value is accepted and pooled into the posterior probability of the parameter value. Repeating these steps many times, we can approximate the probability distribution. A conceptual overview of the ABC rejection sampling algorithm is presented in **Figure 5**.

In the study of colorectal cancer study by Uchi et al. [5], as summary statistics, they adopted the proportions of founder mutations and unique mutations, which is uniquely observed in each sample, in a multiregion mutation profile.

#### **Figure 5.**

*Conceptual overview of the ABC rejection sampling algorithm. This image originally appeared in [18].*

They obtained multiregion mutation profiles for 9 cases with different sample numbers. As the proportions of founder mutations and unique mutations depend on the number of samples, they set the sample number to 5, which is the minimum sample number of the 9 cases, by downsampling the samples in cases containing more than 5 samples. They then estimated the mean of the proportions of founder mutations and unique mutations and used these values as summary statistics values of the observational data (**Figure 6**; note that although we should apply ABC to each of the 9 case separately, they targeted the population mean for simplicity).

For ABC, they generated simulation data while varying 3 parameters, *m* (the mutation rate), *d* (the number of driver genes), and *f* (the effect of driver mutations), which appear to be critical for simulation results (for strategies used to find such parameters, read the next section). In each simulation trial, we simulated multiregion sequencing from a tumor simulated by the BEP model; a multiregion mutation profile was obtained by digging 5 squares out from a simulated tumor and averaging the mutation status of cells in the squares. From the multiregion mutation profile, the proportions of founder mutations and unique mutations were obtained as summery statistics. They performed 50 simulations for each grid point in a threedimensional rectangular parameter space; namely, they assumed a uniform prior for each of the three parameters. For each grid point in the parameter space, they calculate the proportion of the simulation instances whose statistics fall within 1 standard deviation from the mean of the values observed in the real multiregion

*Agent-Based Modeling and Analysis of Cancer Evolution DOI: http://dx.doi.org/10.5772/intechopen.100140*

#### **Figure 6.**

*Fitting the BEP model to multiregion sequencing data by ABC. (A) Observed values of summary statistics in the multiregion sequencing data. After downsampling the samples in cases containing more than 5 samples, the proportions of founder and unique mutations were estimated for 9 cases (case 1–9) and an "average" over the 9 cases was obtained as summary statistics values of the observational data. The error bars at case 1–3 and 5–8 indicate standard deviations for 10 downsampling trials while the error bar at average indicates standard deviations over the 9 cases. (B) Multiregion mutation profiles of 9 colorectal tumors. For the cases except case 4 and 9, representative samples from the downsampling trials were presented as in Figure 4B. This image originally appeared in [5].*

mutation profiles. The distribution of the proportions can be regarded as the posterior and visualized in heat maps (**Figure 7**).

As a result, when cancer evolution was simulated with the assumption of a high mutation rate, we reproduced mutation profiles similar to those obtained by our multiregion sequencing of colorectal cancers (compare **Figure 8A** and **B** with **Figure 4A** and **B**). That is, irrespective of the presence of founder mutations, progressor mutations contributed to the formation of a heterogeneous mutation profile, which was geographically correlated with the sampling locations. Moreover, we also reconstructed local heterogeneity, as illustrated by the finding that progressor mutations existed as mutations with lower allele frequencies in each region. Interestingly, although driver mutations were acquired as founder mutations, progressor mutations contained few driver mutations, and most comprised neutral mutations that did not affect the cell division rate. This suggests that, after the appearance of the common ancestor clone with accumulated driver mutations, extensive ITH was generated by neutral evolution. Moreover, the single-cell mutation profiles of the simulated tumor suggest that the tumor comprises a large number of minute clones with numerous neutral mutations accumulated (**Figure 8C**).

By employing a agent-based model and ABC, Sottoriva et al. [19] also proposed a Big Bang model of human colorectal tumor growth; in their model, tumors grow predominantly as a single expansion producing numerous intermixed subclones that are not subject to stringent selection, which is consistent with the model developed by Uchi *et al.* [5], and both public (clonal) and most detectable private (subclonal) alterations arise early during growth. Hu et al. [20] also employed an agent-based model and ABC to examine the timing of metastasis in colorectal cancer. Multiregion sequencing data containing both primary and metastatic samples were prepared from patients with metastases to the liver or brain. Simultaneously, a spatial agent-based model was developed to simulate tumor growth,

#### **Figure 7.**

*Fitting the BEP model to multiregion sequencing data by ABC (continued). (A) The proportion of simulation instances fitted to the real data. Multiregion mutation profiles were simulated while varying 3 parameters and, for each parameter settings, the proportion of simulation instances that were judged to be similar to the real data based on summary statistics are visualized as heat maps. (B) Multiregion mutation profiles from the simulations. Representative instances from simulation with indicated parameter settings were presented as in Figure 4B. Left blue bars indicate driver genes. This image originally appeared in [5].*

mutation accumulation, and metastatic dissemination. From multiregion sequencing data of each patient, the time of dissemination, which is a parameter in the agent-based model, was estimated by ABC. The results demonstrated that early disseminated cells commonly (81%, 17 of 21 patients) showed metastases, whereas the carcinoma was clinically undetectable (typically, less than 0.01 cm3). Collectively, these examples demonstrated that ABC successfully fitted the simulation model of cancer evolution to cancer genome data, providing insight into the mechanisms of cancer evolution.

Although the problem of computational cost generally accompanies ABC, new sampling approaches utilizing Markov chain Monte Carlo and its derivatives [21]

*Agent-Based Modeling and Analysis of Cancer Evolution DOI: http://dx.doi.org/10.5772/intechopen.100140*

**Figure 8.**

*Computer-simulated tumor with extensive ITH generated by neutral evolution. (A) Tumor simulated based on the BEP model with an assumption of a high mutation rate. (B) Simulated multiregion mutation profile of the simulated tumor. Cell populations in the regions labeled with A–H were extracted from the simulated tumor and their averaged mutation profiles were obtained. (C) Simulated single-cell mutation profile of the simulated tumor. This image was obtained by modifying a figure that originally appeared in [5].*

have been developed to overcome this limitation. Moreover, considering the increasing computing power, this problem will potentially be less important. Notably, ABC has many potential pitfalls [18]. For example, setting the tolerance parameter to zero will give accurate results, but typically at a very high computational cost. In practice, therefore, values of greater than zero are used, but this introduces bias. Similarly, sufficient statistics are sometimes not available and other summary statistics are used instead, but this introduces additional bias because of the loss of information. Additionally, prior distributions and choices of parameter ranges are often subject to criticisms, although they are not unique to ABC and apply to all Bayesian methods. Model complexity (i.e., the number of model parameters) is also an important point. If a model is too simple, it can lack predictive power. In contrast, if the model is too complex, there is a risk of overfitting. Moreover, the complex model faces a problem known as the curse of dimensionality, in which the computational cost is severely increased and may, in the worst case, render the computational analysis intractable. When constructing a simulation model, we should follow the Occam's razor principle: i.e., achieve the lowest model complexity that is sufficient to explain the observational data. To determine the optimal model complexity, we can also employ the model selection scheme based on Bayes factor if a choice of summary statistics is appropriate [22].

#### **4. Characterizing the dynamics of the simulation model**

In the previous section, we explained how to fit a simulation model to observational data. Another direction for studying a simulation is by characterizing the dynamics of the simulation model without observational data. Namely, we can examine parameter dependance by performing a large number of simulations while varying the parameter values. This approach is known as sensitivity analysis and can provide insights into the modeled system as well as identify parameters that are critical for the system dynamics. In sensitivity analysis, as in ABC, we define a summary statistic *Y*. A simulation model is then regarded as a function: *Y* ¼ Fð Þ **X** where **X** ¼ f g *X*1,*X*2, … *Xk* are model parameters. The aim of sensitivity analysis can also be considered as characterizing the function "F".

So far, a number of approaches have been proposed for sensitivity analysis. For example, one-factor-at-a-time (OFAT) sensitivity analysis is one of the simplest and most common approaches that changes one parameter at a time to determine

the effects on a summery statistic [23]. In OFAT sensitivity analysis, we move one parameter, while leaving the other parameters at their baseline (nominal) values, and then return the parameter to its nominal, which is repeated for each of the other parameters. We then plot the relationship between each parameter and a summary statistic to examine the dependency of the summary statistic on the parameter, or the relationship can be measured by partial derivatives or linear regression. In exchange for its simplicity, this approach does not fully explore the input space, as it does not consider the simultaneous variation of multiple parameters. This means that the OFAT approach cannot detect interactions between parameters.

Global sensitivity analysis aims to address this point by sampling a summary statistic over a wide parameter space involving multiple parameters. Sobol's method is a popular approach for estimating the contributions of different combinations of parameters to the variance of the summary statistic while assuming that all parameters are independent [24]. The sensitivity of the summary statistic *Y* to a parameter *Xi* is measured by the amount of variance in *Y* caused by the parameter *Xi* and can be expressed as a conditional expectation, Var Eð Þ **X**�*i* ð Þ *Y*j*Xi* , where "V*ar*" and "E" denote the variance and expected value operators, respectively, and **X**�*i* denotes the set of all input variables except for *Xi*. This expression essentially measures the contribution *Xi* alone to uncertainty (variance) in *Y* (averaged over variations in other variables), and is known as the first-order sensitivity index or main effect index. Importantly, it does not measure the uncertainty caused by interactions with other variables. A further measure, known as the total effect index, gives the total variance in *Y* caused by *Xi* and its interactions with any of the other input variables. Both quantities are typically standardized by dividing by Varð Þ *Y* . In Sobol's method, we typically attempt full exploration of the parameter space based on a Monte Carlo method to grasp parameter interactions and nonlinear responses.

However, such approaches appears to be insufficient to comprehensively grasp how the parameters judged to be influential control the behaviors of agent-based models. To overcome this point, Niida *et al.* [25] recently developed a new approach to sensitivity analysis for agent-based simulations, named as MASSIVE (Massively parallel Agent-based Simulations and Subsequent Interactive Visualization-based Exploration). MASSIVE overcomes the limitations of existing methods by taking advantage of two currently rising technologies: massively parallel computation and interactive data visualization. MASSIVE employs a full factorial design involving a multiple number of parameters (i.e., test every combination of candidate values of the multiple parameters), which can broadly cover a target parameter space. In addition, when analyzing a stochastic simulation model such as an agent-based model, multiple simulation trials with the same parameter setting are required to examine stochastic effects. To cope with the computational cost problem caused by these features, MASSIVE utilizes a supercomputer, in which agent-based simulations with different parameter settings and the following post-processing step of simulation results are performed in parallel. The massively parallel simulations generate massive results, which then pose a problem for interpretation. This problem was solved by developing a web-based tool that interactively visualizes not only the values of multiple summary statistics, but also output images (e.g., mutation profiles) from simulations with each parameter setting.

Below I explain an example of sensitivity analysis, which was performed by Niida et al. [26] to understand the precise mechanisms underlying neutral evolution induced by a high mutation rate. First, they built an agent-based model, referred to as the "neutral" model, for simulating neutral evolution in cancer. Although the neutral model is similar to the BEP model, the neutral model assumes only neutral mutations and omits spatial information. They also improved the approach used for mutation accumulation in the BEP model. Namely, in the neutral model, they

considered only neutral mutations that did not affect cell division and death. In a unit time, a cell divides into two daughter cells with a constant probability *g*0 without dying. In each cell division, each of the two daughter cells acquires *kn* � Poisð Þ *mn=*2 neutral mutations. They assumed that neutral mutations acquired by different division events occur at different genomic positions.The simulation started from one cell without mutations and ended when the population size *p* reached *P* or time *t* reached *T*.

Through sensitivity analysis based on the MASSIVE method, they confirmed that the mutation rate is the most important factor affecting neutral evolution (**Figure 9**). As a summary statistic for evaluating ITH, they calculated Shannon index 0.05 by the following procedure. After removing mutations with frequencies of less than 0.05, the proportions of different subclones (cell subpopulations with different mutations) were obtained and the Shannon index *H* was calculated using the following formula: *H* ¼ �P*n i*¼1*pi* log *pi* � �, where *n* is the total number of different subclones and *pi* is the proportion of each subclone. Based on this definition, a larger Shannon index 0.05 value indicates more extensive ITH. Together with a heat map of the Shannon index 0.05 values, we visualized single-cell mutations profiles obtained for different parameter settings. The mutation profile matrix was obtained by sampling 1,000 cells from a simulated tumor, and visualized after filtering out lower-frequency mutations, such that the maximum number of rows was 300. The rows and columns are reordered by hierarchical clustering and index mutations and samples, respectively. They found that when the mean number of mutations generated by per cell division, *mn*, was less than 1, the neutral model just generated sparse mutation profiles with relatively small values of Shannon index 0.05. In contrast, when *mn* exceeded 1, the mutation profiles presented extensive ITH, which are characterized by a fractal-like pattern and large values of the ITH score (hereinafter, this type of ITH is referred to as "neutral ITH"). These results suggest that neutral ITH is shaped by neutral mutations that trace the cell lineages in the simulated tumors. Note that the mutation profiles were visualized after filtering out low-frequency mutations. Assuming a high mutation rate, more numerous subclones with different mutations should be observed if mutations existing at lower frequencies are counted. However, the ITH score does not depend on the population size *P* because low-frequency mutations were filtered out before calculation.

Thus far, several theoretical and computational studies have shown that a stem cell hierarchy can boost neutral evolution in a population of cancer cells [12, 27];

#### **Figure 9.**

*Sensitivity analysis of the neutral model. (A) Heap map obtained by calculating* Shannon index 0.05 *while changing the neutral mutation rate mn and maximum population size P. (B–H) Single-cell mutations profiles obtained for seven parameter settings, which are indicated on the heat map in A. This image originally appeared in [26].*

based on this, they extended the neutral model to the "neutral-s" model such that it contains a stem cell hierarchy (**Figure 10**). The neutral-s model assumes that two types of cell exist: stem and differentiated. Stem cells divide with a probability *g*0 without dying. For each cell division of stem cells, a symmetrical division generating two stem cells occurs with a probability *s*, whereas an asymmetrical division generating one stem cell and one differentiated cell occurs with a probability 1 � *s*. A differentiated cell symmetrically divides to generate two differentiated cells with a probability *g*0 but dies with a probability *dd* 0. The means of accumulating neutral mutations in the two types of cell is the same as that in the original neutral model, which means that the neutral-s model is equal to the original neutral model when *s* ¼ 0 or *dd* 0 ¼ 0. For convenience, they define *δ* ¼ log 10 *dd* 0*=g*0 and hereinafter use *δ* rather than *dd* 0.

MASSIVE analysis of the neutral-s model confirmed that incorporation of the stem cell hierarchy boosts neutral evolution (**Figure 11**). To obtain the heat map in **Figure 11A**, the ITH score was measured while *dd* 0 and *δ* were changed, whereas *mn* ¼ 0*:*1 and *P* ¼ 1000 were maintained as constant. In the heat map, a decrease in *s* leads to an increase in the ITH score when *δ*≥ 0 (i.e., *dd* 0 ≥*g*0). A smaller value of *s* means that more differentiated cells are generated per stem cell division, and *δ*≥0 means that the population of differentiated cells cannot grow in total, which is a valid assumption for typical stem cell hierarchy models. That is, this observation indicates that the stem cell hierarchy can induce neutral ITH even with a relatively low mutation rate setting (i.e., *mn* ¼ 0*:*1), with which the original neutral model cannot generate neutral ITH.

The underlying mechanism boosting neutral evolution can be explained as follows. Only stem cells were considered for an approximation, as differentiated cells do not contribute to tumor growth with *δ*≥0. While one cell grows to a population of *P* cells, let cell divisions synchronously occur across *x* generations during the clonal expansion. Then, 1ð Þ þ *s x* ¼ *P* holds because the mean number of stem cells generated per cell division is estimated as 1 þ *s*. Solving the equation for *x* gives *x* ¼ log *P=* log 1ð Þ þ *s* ; that is, it can be estimated that during clonal expansion, each of the *P* cells experiences log *P=* log 1ð Þ þ *s* cell divisions and accumulates *mn* log *P=*2 log 1ð Þ þ *s* mutations on average. They confirmed that the expected mutation count based on this formula fit well with the values observed in their simulation (data not shown). These arguments mean that a tumor with a stem cell hierarchy accumulates more mutations until reaching a fixed population size than

#### **Figure 10.**

*Schema of the neutral-s model. Stem cells divide with a probability go without dying. For each cell division of stem cells, a symmetrical division generating two stem cells occurs with probability s, whereas an asymmetrical division generating one stem cell and one differentiated cell occurs with probability* 1 � *s. A differentiated cell symmetrically divides to generate two differentiated cells with probability g*0 *but dies with probability dd* 0*. This image originally appeared in [26].*

*Agent-Based Modeling and Analysis of Cancer Evolution DOI: http://dx.doi.org/10.5772/intechopen.100140*

#### **Figure 11.**

*Sensitivity analysis of the neutral-s model. (A) Heat map obtained by calculating Shannon index 0.05 while changing the relative death rate of differentiated cells δ* ¼ log 10 *dd* 0*=g*0 *and symmetrical division rate s. The neutral mutation rate mn and maximum population size P set to* 10�1 *and* 105*, respectively. (B–J) Single-cell mutation profiles obtained for nine parameter settings, indicated on the heat map presented in A. This image originally appeared in [26].*

does a tumor without a stem cell hierarchy. That is, a stem cell hierarchy increases the apparent mutation rate by log 2*=* log 1ð Þ þ *s* -fold, which induces neutral evolution even with relatively low mutation rate settings.

Recent genomic analysis demonstrated that multiple evolutionary modes exists in cancer systems. For example, as described above, ITH in renal cancer is generated by Darwinian selection, which is in contrast to neutral evolution in colorectal cancer. Moreover, by multiregion sequencing of early-stage colorectal tumors, Saito *et al.* [26] showed that ITH is shaped by Darwinian selection in the early phase of colorectal cancer evolution, which means that a temporal shift of the evolutionary principle shaping ITH occurs during colorectal tumorigenesis. Employing agentbased modeling and MASSIVE analysis, Niida *et al.* [26] also constructed a model that explain this evolutionary shift. Darwinian ITH in an early-stage tumor is reproduced by the assumption of multiple driver mutations of relatively weak strength. At some point, growth of the early colorectal tumor slows because resource limitations, which is reproduced by introducing the carrying capacity into the simulation model. When they assumed that an explosive mutation that negates the carrying capacity was obtained with a small probability, a clone acquiring the explosive mutation overcame the resource limitation and expanded as late-stage tumors, in which ITH was generated neutral evolution. Another simulation study by West *et al.* [28] proposed that spatial constraints and limited cellular mixing play important roles in a similar Darwinian-neutral shift.

Sensitivity analysis also provides insight into metastatic tumor progression, which is poorly understood despite its clinical importance. Evaluation of genomic divergence between paired metastatic and primary tumors (M-P divergence) from multiregion sequencing is a good starting point for addressing this problem. Sun and Nikolakopoulos [29] extended *tumopp* [14] to simulate paired primary and metastatic tumors, and explored factors affecting M-P divergence by sensitivity analysis. As a result, they found that M-P divergence depends not only on the metastatic dissemination time, but also on the evolutionary dynamics and

detectability of seedling cell lineages in a primary tumor. It was concluded that investigating tumor growth dynamics in detail is important, particularly when researchers interpret heterogeneity among longitudinal samples to infer the evolutionary timeline of cancer progression. Collectively, these examples demonstrated that agent-based modeling combined with sensitivity analysis is a useful tool for studying cancer evolutionary dynamics.


1. Introduction

There is an ongoing in the theory of interaction of charged particle beams with plasmas. Although most theoretical works have reported on the energy loss of ions in a plasma without magnetic field, the strongly magnetized case has not yet received as much attention as the field-free case. The energy loss of ion beams and the related processes in magnetized plasmas are important in many areas of physics such as transport, heating, magnetic confinement of thermonuclear plasmas, and astrophysics. The range of the related topics includes ultracold plasmas [1, 2], the

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited. © 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

cooling of heavy ion beams by electrons [3–12], as well as many very dense systems involved in magnetized target fusions [11], or heavy ion inertial confinement fusion (ICF).

with corresponding Larmor radius indicates that to experience a strong influence of the

Stopping Power of Ions in a Magnetized Plasma: Binary Collision Formulation

endeavors to lead to the very first unambiguous and genuine identification of an experimental

Motivated by these recent developments, our purpose is to investigate the SP of an ion moving in a magnetized plasma in a wide range of the value of a steady magnetic field. The present paper is based on our earlier studies in Refs. [8, 24, 44, 45] where the second-order energy transfers for individual collisions of electron-ion [8, 24, 44] of any two identical particles, like electron-electron [44], and finally of two gyrating arbitrary charged particles [45] have been calculated with the help of an improved BC treatment. This treatment is—unlike earlier approaches of, e.g., Refs. [9, 42]—valid for any strength of the magnetic field. As the first application of the theoretical BC model developed in Refs. [8, 24, 44, 45], we have calculated in Ref. [47] the cooling forces on the heavy ion beam interacting with a strongly magnetized and temperature anisotropic electron beam. It has been shown that there is a quite good overall agreement with both the CTMC numerical simulations and the experiments performed at the ESR storage ring at GSI [48–50].

In Section 2 we introduce briefly a perturbative binary collision formulation in terms of the binary force acting between an ion and a magnetized electron and derive general expressions for the second-order (with respect to the interaction potential) stopping power. In contrast to the previous investigations in Refs. [8, 24, 44, 45], we here consider the (macroscopic) stopping force which is obtained by integrating the binary force of an individual electron-ion interaction with respect to the impact parameter and the velocity distribution function of electrons. That is, the stopping force for monoenergetic electrons is folded with a velocity distribution. The resulting expressions involve all cyclotron harmonics of the electrons' helical motion and are valid for any interaction potential and any strength of the magnetic field. In Section 2.4 we present explicit analytic expressions of this second-order stopping power for the specific case of a regularized and screened interaction potential [51, 52] which is both of finite range and less singular than the Coulomb interaction at the origin and which includes as limiting cases the Debye (i.e., screened) and the Coulomb potentials. For comparison of our expressions with previous approaches, we consider in Section 3 the corresponding asymptotic expressions for large and small ion velocities and strong and vanishing magnetic fields. The analytical expressions presented in Section 2.4 are evaluated numerically in Section 4 using parameters of the envisaged experiments on ion stopping [46]. In particular, we compare our approach with the CTMC simulations. The results are summarized and discussed in Section 5. The regularization parameter and the screening length involved in the interaction potential are briefly specified and discussed in Appendix A.

Let us consider two point charges with masses m, M and charges �e, Ze, respectively, moving in a homogeneous magnetic field B ¼ Bb. We assume that the particles interact with the potential �Z=e 2U rð Þ with =e 2 ¼ e2=4πe0, where e0 is the permittivity of the vacuum and

. We expect these

69

http://dx.doi.org/10.5772/intechopen.77213

magnetic field, the electron density should be comparable with a few 1016 cm3

magnetic signature for nonrelativistic ion stopping in plasmas.

2. Theoretical model

2.1. Binary collision (BC) formulation

For a theoretical description of the energy loss of ions in a plasma, there exist some standard approaches. The dielectric linear response (LR) treatment considers the ion as a perturbation of the target plasma, and the stopping is caused by the polarization of the surrounding medium. It is generally valid if the ion couples weakly to the target. Since the early 1960s, a number of calculations of the stopping power (SP) within LR treatment in a magnetized plasma have been presented (see Refs. [13–37] and references therein). Alternatively, the stopping is calculated as a result of the energy transfers in successive binary collisions (BCs) between the ion and the electrons [37–45]. Here, it is necessary to consider appropriate approximations for the screening of the Coulomb potential by the plasma [8]. However, significant gaps between these approaches involve the ion stopping along magnetic field B and perpendicular to it. In particular, at high B values, the BC predicts a vanishingly parallel energy loss, which remains at variance with the nonzero LR one. Also, challenging BCLR discrepancies persist in the transverse direction, especially for vanishingly small ion projectile velocity vi when the friction coefficient contains an anomalous term diverging logarithmically at vi ¼ 0 [23, 24]. For calculation of the energy loss of an ion, two new alternative approaches have been recently suggested. One of these methods is specifically aimed at a low-velocity energy loss, which is expressed in terms of velocity-velocity correlation and, hence, to a diffusion coefficient [34]. Next, in Ref. [27] using the Bhatnagar-Gross-Krook approach based on the Boltzmann-Poisson equations for a collisional and magnetized classical plasma, the energy loss of an ion is studied through a LR approach, which is constructed such that it conserves particle number locally.

An alternative approach, particularly in the absence of any relevant experimental data, is to test various theoretical methods against comprehensive numerical simulations. This can be achieved by a particle-in-cell (PIC) simulation of the underlying nonlinear Vlasov-Poisson Equation [10, 31]. While the LR requires cutoffs to exclude hard collisions of close particles, the collectivity of the excitation can be taken into account in both LR and PIC approaches. In the complementary BC treatment, the stopping force has been calculated numerically by scattering statistical ensembles of magnetized electrons from the ions in the classical trajectory Monte Carlo (CTMC) method [7, 10, 37–41]. For a review we refer to a recent monograph [8] which summarizes all theoretical and numerical methods and approaches also discussing the ranges of their validity.

The very recent upheaval of successful experiments involving hot and dense plasmas in the presence of kilotesla magnetic fields (e.g., at ILE (Osaka), CELIA (Bordeaux), LULI (Palaiseau), LLNL (Livermore)) remaining nearly steady during 10–15 ns strongly motivates the fusion as well as the warm dense matter (WDM) communities to investigate adequate diagnostics for their dynamic properties. This opens indeed a novel perspective by allowing magnetic fields to play a much larger if not a central role both in ICF and WDM plasmas. In this context proton or any nonrelativistic ion stopping is likely to provide an option of choice for investigating genuine magnetization features such as anisotropy, when the electron plasma frequency turns significantly lower than the cyclotron one [46]. In addition, an experimental test of proton or alpha particle stopping in a magnetized plasma is currently envisioned (see, e.g., Ref. [46] for a preliminary discussion). The parameters at hand are a fully ionized hydrogen plasma with a density up to 1020 cm3 and temperature between 1 and 100 eV. The steady magnetic field can be up to 45 T strong. A preliminary examination based on comparing electron Debye length with corresponding Larmor radius indicates that to experience a strong influence of the magnetic field, the electron density should be comparable with a few 1016 cm3 . We expect these endeavors to lead to the very first unambiguous and genuine identification of an experimental magnetic signature for nonrelativistic ion stopping in plasmas.

Motivated by these recent developments, our purpose is to investigate the SP of an ion moving in a magnetized plasma in a wide range of the value of a steady magnetic field. The present paper is based on our earlier studies in Refs. [8, 24, 44, 45] where the second-order energy transfers for individual collisions of electron-ion [8, 24, 44] of any two identical particles, like electron-electron [44], and finally of two gyrating arbitrary charged particles [45] have been calculated with the help of an improved BC treatment. This treatment is—unlike earlier approaches of, e.g., Refs. [9, 42]—valid for any strength of the magnetic field. As the first application of the theoretical BC model developed in Refs. [8, 24, 44, 45], we have calculated in Ref. [47] the cooling forces on the heavy ion beam interacting with a strongly magnetized and temperature anisotropic electron beam. It has been shown that there is a quite good overall agreement with both the CTMC numerical simulations and the experiments performed at the ESR storage ring at GSI [48–50].

In Section 2 we introduce briefly a perturbative binary collision formulation in terms of the binary force acting between an ion and a magnetized electron and derive general expressions for the second-order (with respect to the interaction potential) stopping power. In contrast to the previous investigations in Refs. [8, 24, 44, 45], we here consider the (macroscopic) stopping force which is obtained by integrating the binary force of an individual electron-ion interaction with respect to the impact parameter and the velocity distribution function of electrons. That is, the stopping force for monoenergetic electrons is folded with a velocity distribution. The resulting expressions involve all cyclotron harmonics of the electrons' helical motion and are valid for any interaction potential and any strength of the magnetic field. In Section 2.4 we present explicit analytic expressions of this second-order stopping power for the specific case of a regularized and screened interaction potential [51, 52] which is both of finite range and less singular than the Coulomb interaction at the origin and which includes as limiting cases the Debye (i.e., screened) and the Coulomb potentials. For comparison of our expressions with previous approaches, we consider in Section 3 the corresponding asymptotic expressions for large and small ion velocities and strong and vanishing magnetic fields. The analytical expressions presented in Section 2.4 are evaluated numerically in Section 4 using parameters of the envisaged experiments on ion stopping [46]. In particular, we compare our approach with the CTMC simulations. The results are summarized and discussed in Section 5. The regularization parameter and the screening length involved in the interaction potential are briefly specified and discussed in Appendix A.


As visible in Fig. 5.19, the set X o 20 leads to the fastest approach. Since the reactive decelerations along with their probabilities are considered comfortable, this approach is pursued. The considered decelerations include those to belief updates w.r.t. certainly existing objects and also those to potentially occurring objects behind the occlusion. The runtime of the planner for this scenario with the presented belief progressions on a single core of an Intel® Core™ i7-8565U was below 100 ms. Parallelized computation for the different belief progressions could further lower the computation time. Yet, runtime analysis or optimization is not the focus of this work. To summarize, the approach works well in presence of many objects, given that the prediction module provides suitable belief progressions.

(a) Attempted velocity profile, assuming obj3 goes straight and obj2 does not exist (X o 20). The reactions to the possible belief updates are depicted dashed. The two reactions that lead to a full stop in front of the CZ (X o 21,X o 22) are almost alike, as they start only 1 s apart from each other (slight deviations in a from 2 s to 6 s).

(b) Alternative velocity profile, assuming obj3 turns right (X o 10,egoalt10). Since certainty about this right turn is only reached once it turns, a change to going straight is not neglected in the safety checks, leading to a slow down compared to the attempted profile.

(c) Alternative velocity profile, assuming obj3 turns left (X o 30,egoalt30). The ego vehicle has to slow down further to merge in behind obj3.

Figure 5.19: Attempted velocity profile and alternatives for the scenario from Fig. 5.18. Since sets are considered, the objects' positions are depicted with an area instead of a line.

### **5.2.2 Lane Change**

In this section, the approach to courtesy-aware lane changes from Sec. 4.4 is evaluated in an urban lane change scenario.

#### **Implementation Details**

The implementation presented in [9] was integrated into a lanelet2-based routing state machine and the previously explained modified control (cf. Sec. 5.1) has been employed. Further, cases with only one or no vehicle in the target lane have been added to the lane change state machine. As heuristic for the gap selection, the gap that can be reached earliest in space along the target lane with a trajectory within the comfortable acceleration limits is selected. The approach to the gap is computed analytically as a sequence of constant acceleration, constant velocity and constant acceleration. If the gap is currently ahead, regarding position and velocity difference, the vehicle accelerates up to a certain velocity v + peak to catch up, before decelerating to the gap velocity. If the gap is behind, the vehicle decelerates to a velocity v − peak before accelerating to the gap velocity. The path for the lane change is computed as transition from the path in the source lane to the path in the target lane using the sigmoid function tanh. The planning parameters can be found in Appendix A.2.

#### **Evaluation Scenario**

The general performance of the approach to facilitate lane changes is shown in a mandatory lane change at Haid-und-Neu-Straße in Karlsruhe, Germany (cf. Fig. 5.20), as presented in [9]. The success of approaches relying on cooperation or courtesy of other traffic participants inherently depends on whether the respective agents are willing to cooperate or not. Instead of providing velocity or acceleration profiles for this multi-lane scenario, we perform a qualitative evaluation of the lane change behavior: If the vehicle behind us in the target lane is courteous, the lane change is assumed to be safe after the additional lane change response time has been granted and thus can be performed. Courtesy of the vehicle behind was modeled using the enhanced IDM [KTH10] for the motion planning of this vehicle and considering the ego vehicle as front vehicle as soon as it indicated its lane change desire. In case the respective vehicle does not open a gap of sufficient size for RSS safety, even though the desire to change the lane has been communicated using the indicator and moving left within the current lane, the lane change is aborted after waiting for ∆*t*wait. The non-courteous behavior was modeled by refraining from considering vehicles in other lanes as front vehicle.

This experiment shows the general applicability to scenarios where approaches with an independent prediction would lead to overly conservative behavior, i.e. waiting for large gaps without "asking for them" via clear communication of the lane change desire. Besides this advantage in comfort and convenience, safety is ensured w.r.t. the defined assumptions of the presented RSS concept for *non-reckless* lane changes. Since the reaction time for the ego vehicle can be significantly lower than the commonly granted one second for human drivers, the safety distance to the front vehicle can be significantly smaller than the one that has to be granted to the vehicle behind (cf. Fig. 5.21). This leads to less additional space needed for lane changes. Both lane change attempts are also shown in the video5 provided along with [9].

The runtime of the planner for this scenario on a single core of an Intel® Core™ i7-8565U was below 2 ms for all plans throughout the simulation. While runtime analysis or optimization is not the focus of this work, the computation time underlines the general suitability of the approach for realtime planning in automated vehicles. In order to further show the real-time capability and applicability of the approach, it has also been implemented on our experimental vehicle *BerthaOne*.

5 Also available via http://www.mrt.kit.edu/z/publ/download/2019\_LaneChange\_Naumann.mp4, last retrieved 2020-05-28.

Figure 5.20: Aerial image of the lane change scenario. Imagery © 2020 Google, map data © 2020 Stadt Karlsruhe VLW, GeoBasis-DE/BKG (©2009).

Figure 5.21: Screenshot of the lane change in the simulation framework CoInCar-Sim [5] (from [9], © 2019 IEEE). The other vehicles are depicted in blue. The ego vehicle is depicted in black. Its trajectory is visualized with the colored circles at equidistant times.

### **5.2.3 Mutual Cooperation**

In addition to the sampling-based implementation presented in [3, 6], we show the performance of the approach to mutual cooperation from Sec. 4.5 with a C++ implementation in a narrowing scenario.

#### **Implementation Details**

In contrast to the implementation presented in [3, 6], the global optimum is no longer found via sampling, but via employing the enhanced IDM [KTH10]. For the globally optimal trajectory ensemble, vehicles are no longer expected to accelerate beyond the velocity that they would have chosen if no other vehicle was present. The actual human behavior in cooperative situations is to be investigated further by experts in this field and is not part of this work. Yet, this second implementation shows the flexibility of the approach regarding different models of cooperation.

In order to investigate whether cooperation is possible, i.e. whether the globally optimal solution is not ambiguous, the homotopy classes of the solutions for different driver models are analyzed. In case they unambiguously yield that the ego vehicle should go first, a predefined model compliance probability *p*mc is assumed. As model compliance check for other drivers, the deviation from the expected velocity ∆vmc is considered. A model compliance violation is expected to be measurable by the time when the predictions for the different homotopy classes of the cooperative vehicle differ by ∆vmc. For this point in time, we compute the necessary deceleration to switch to going second. With this information, the heuristic presented in Sec. 4.3 is employed to check whether attempting the globally optimal solution violates the comfort margin, as the expected deceleration for the fallback plan is too high. In case the comfort margin is not violated, the globally optimal solution is attempted. The expected velocity profile for the other vehicle is stored.

If the ego vehicle is still far away from the conflict zone, however, the situation might be currently ambiguous, or the velocity difference in-between the homotopy classes might be not detectable, but the situation might become clear later. Thus, it is checked whether continuing in the hope of the cooperative solution is convenient. Since the time and the probability of a positive decision for the cooperative solution are hard to estimate, the decision is only postponed if (a) the necessary deceleration for going second in the next step is below a certain value *a*decel,low or (b) the maximum deceleration remains unchanged, for example as a deceleration due to significant curvature is upcoming anyway.

If neither attempting to go first, nor decision postponing is considered comfortable, the ego vehicle attempts to go second. In case the model compliance is violated, i.e. the velocity of the cooperative vehicle is more than ∆vmc higher than expected, the ego vehicle switches to go second as the globally optimal solution is no longer believed to be realized.

The planning parameters can be found in Appendix A.2.

#### **Evaluation Scenario**

The approach is evaluated at an artificial narrowing in Karlsruhe (cf. Fig. 5.22, Fig. 5.23). The vehicle configuration is chosen such that the desired trajectories for both vehicles, each assuming the other vehicle was not there, would result in a collision inside the conflict zone (cf. Fig. 5.24). For such scenarios, mutually predicting the other vehicle independently, for example with constant velocity or with the IDM, commonly results in a deadlock, at least until one vehicle comes to a full stop. Fig. 5.25 shows the planning result of both vehicles with independent prediction after 8 s. A similar approach was for example used by Ziegler et al. [ZBS+14]. With the presented approach, cooperation is possible in this scenario, in case the other driver behaves cooperatively. Meaningful malicious behavior, i.e. egoistic behavior, can also be detected and dealt with comfortably (cf. Fig. 5.26). Adversarial behavior of others, such as provoking deadlocks, is still possible, but not expected and not the focus of this work. The runtime of the planner for this scenario on a single core of an Intel® Core™ i7-8565U was below 20 ms for all plans throughout the simulation. While runtime analysis or optimization is not the focus of this work, the computation time underlines the general suitability of the approach for real-time planning in automated vehicles.

Figure 5.22: Aerial image of the narrowing scenario. Imagery © 2020 Google, map data © 2020 Stadt Karlsruhe VLW,GeoBasis-DE/BKG (©2009).

Figure 5.23: Paths in the narrowing scenario.

Figure 5.24: Possible solutions from the ego perspective: Expected is the solution depicted in solid, with the ego vehicle going first. The trajectory for the other vehicle acting egoistic is depicted in brown dashed. The reaction to this behavior, by the time a model violation is expected to be detectable, is depicted in blue dashed. The alternative of assuming straightaway that the other will act egoistic, and thus attempting to go second, is depicted in light blue dashed.

Figure 5.25: Trajectories for mutually predicting each other with constant velocity: The already driven part is visualized with a solid line, the prediction for the other vehicle is dotted, the own plan is dashed.

Figure 5.26: Driven trajectories with the proposed approach: In case obj1 behaves cooperatively as expected, the ego vehicle drives first (solid). In case obj1 acts egoistically, the ego vehicle detects this and reacts early and comfortably (dashed).

# **6 Conclusions and Future Directions**

# **6.1 Conclusions**

In this work, several contributions to motion planning for automated vehicles in mixed traffic have been presented. At first, the different sources of uncertainties and common problem formulations in the context of motion planning for automated vehicles are outlined. This consideration exposed that, while optimal control approaches with one deterministic trajectory as result are applicable, they can be outperformed by approaches that account for the repeated replanning. Particularly in cases where the possible future evolutions of a scene largely differ, for example due to multiple route options for other traffic participants, considering multiple future trajectories is beneficial. Yet, probabilistic approaches that attempt a comprehensive uncertainty treatment with limited runtime are often operating on a largely restricted action space. Also, many approaches do not consider safety and compliance with the traffic rules.

Regarding the safety verification, recent frameworks focus on the notion of responsibility. That is, instead of guaranteeing to never be involved in an accident, automated vehicles should guarantee that they will never cause an accident. The traffic rules and the right of way defined therein have to be accounted for when defining responsibility. Thus, the safety of a motion plan does not only depend on the road layout and the current vehicle configurations, but also on the priorities of lanes and roads over each other. This motivates the decomposition of the safety concept into scenarios: As a contribution to safety in motion planning, in Sec. 4.2, the existing responsibility sensitive safety (RSS) framework [SSS18] is enhanced by the consideration of traffic rules at intersections with crossing and merging traffic, including occlusions and limited sensor range, as well as lane changes. In particular, the consideration of right-of-way rules at intersections leads to an upper bound for the velocity to ensure safety in presence of occlusions and limited sensor range.

Given these hard constraints from safety and the traffic rules, the goal is to generate convenient motion plans, i.e. plans that are comfortable but not overcautious. Towards this goal, this work relies on the assumption that, due to the complexity of road traffic, different scenarios can comprise inherently different types of uncertainties that have to be taken into account. These differences largely affect the suitability of different problem formulations. While some uncertainties have to be actively reduced, others can hardly be influenced, for example. Thus, the main contribution of this work is the introduction of three motion planning approaches, targeted towards the different predominant uncertainties in different scenarios.

For scenarios with unimpeded sight and unambiguous routes for all traffic participants, where mutual cooperation is required, a global optimum approach with model compliance consideration is presented in Sec. 4.4. The approach is based on a multi-agent formulation, where other vehicles are expected to cooperate and to have a similar understanding of a globally optimal solution. At the same time, it is checked whether the global optimum is unique, also for different measures of optimality, being aware that the cooperation model might be violated. In case the global optimum is unambiguous and the other traffic participant behaves within our model bounds, cooperation is facilitated, resulting in globally optimal trajectories. Otherwise, the ego vehicle reacts early and pursues a comfortable fallback solution.

For lane change scenarios, in which courtesy of others is required in case of dense traffic, a second approach is presented in Sec. 4.5. The approach is rule-based, and also explicitly modeling other drivers' behavior as part of the decision problem. Even for gaps that are currently too narrow for a safe lane change, the lane change desire is indicated and the response of the potentially courteous vehicle is awaited. In case this vehicle opens a gap for us that is large enough for a safe lane change according to the enhanced RSS framework, the lane change is pursued. Otherwise, the next available gap is to be approached. The approach facilitates the design of gap approach and desire indication according to results from human behavior studies. Furthermore, to the best of the author's knowledge, the presented approach is the first lane change approach that examines safety also w.r.t. unexpected but lawful human behavior, such as emergency decelerations or closing a gap that would have allowed a lane change. It facilitates lane changes on urban multi-lane roads as well as on highways.

Large parts of the urban and rural road network consist of a single lane per driving direction and intersections with predefined right of way. The most significant contributions of this thesis are the formulation of this problem as a partially observable Markov decision process (POMDP) with an ego-actionindependent belief over the trajectories of other traffic participants in Chapter 3 and the presentation of an approach to simplify and solve this problem formulation in Sec. 4.3. The approach generates trajectories based on uncertainty estimates for the trajectories of other traffic participants. The uncertainties include their current state and their desired route, for example. Not only the current uncertainty estimate, but also possible future progressions are incorporated. Furthermore, occlusions and limited sensor range, along with an estimate of their future evolution, and even potential ghost objects or estimates of objects behind the sensing horizon are considered. Also, appropriate reactions to uncertain traffic light state classifications or pedestrians that walk towards a crosswalk can be determined by the approach. Safety is ensured based on the enhanced RSS framework. Related approaches that solve this POMDP in a belief space with algorithms based on Monte Carlo tree search (MCTS) build an action-dependent belief tree during planning. In this process, knowledge about the behavior of others is commonly modeled via a probabilistic transition model, hidden states, and a probabilistic observation model. The assumption behind the presented approach is that the behavior of prioritized vehicles is independent of the ego behavior. The observations that facilitate a behavior prediction are only influenced by occlusions and the sensor range, which again depend on the position of the ego vehicle along its route. Yet, the belief tree can be built independently of the ego motion plan in an upstream prediction module. Here, more sophisticated behavior and observation models can be employed, since they are not repeatedly queried during planning. Instead, only the result of those sophisticated models, provided through selected belief progressions, is used during planning. To the best of the author's knowledge, the presented approach is the first to comprehensively consider uncertainties regarding the current state and future trajectory of other traffic participants, including multiple routes, and regarding their existence, based on object detections within the field of view, but also beyond the sensor range and in occluded areas, while operating in a continuous action space and while guaranteeing safety w.r.t. to the enhanced RSS framework. Lastly, the approach is suitable for scenarios with several traffic participants, since the combinatorial burden can be largely reduced through the limitation to selected belief progressions.

# **6.2 Future Directions**

On the way towards safe and convenient motion planning for automated vehicles in mixed traffic, several research directions are promising. The presented approaches can be extended to scenarios that are subject to special traffic laws, such as the zipper merge, or overtaking with oncoming traffic.

Further, to allow for more efficient cooperation, the models of human behavior can be improved using naturalistic driving datasets, such as the recently presented INTERACTION dataset [10]. Fundamental work on how to reveal cost functions from parts of human-driven trajectories has been presented in [11].

Another research direction is a planning-focused scene understanding and prediction. Computing a full trajectory for all potential trajectories of others, along with the fallback trajectories in case of deviations, is infeasible. Thus, the behavior of others was grouped into belief progressions with a set-based motion prediction. Here, the right granularity has to be found as a trade-off between keeping the problem formulation feasible and solvable in real-time, while not discarding promising maneuver ensembles. With this input, the presented belief-dependent planning approach can be evaluated in experimental vehicles.

Further work in the area of scene understanding and prediction includes heuristics for cooperation probabilities with other traffic participants, in order to improve the gap choice for desired lane changes, but also to get a better confidence estimate for mutually cooperative scenarios.

# **A Appendix**

## **A.1 Safe Distance for Merging**

As stated in eq. 4.5, for the time *t* > 0 after the front vehicle has merged, the safety reserve for the prioritized rear vehicle to the front vehicle is

$$d\_{\rm reservoir}(t) = s\_{\rm merge}^{-}(t) - s\_{\rm prior}^{+}(t) - d\_{\rm safe,prio}(t). \tag{A.1}$$

It is assumed that the rear vehicle no longer accelerates from the time the front vehicle merges, visualized in Fig. 4.2, as this maneuver is foreseeable. Still, the usual reaction time ρprio is granted before the vehicle is expected to start decelerating.

For *t* > ρprio, the position of the front bumper of the prioritized vehicle is

$$s\_{\rm proj}^{+}(t) = -d\_{\rm init,prio}^{+} + \nu\_{0,prio}t + \frac{1}{2}a\_{0,prio}(t - \rho\_{\rm proj})^2. \tag{A.2}$$

The position of the rear bumper of the merging vehicle is

$$s\_{\rm merge}^{-}(t) = -d\_{\rm init,\rm merge}^{-} + \nu\_{0,\rm merge}t + \frac{1}{2}a\_{0,\rm merge}t^2. \tag{A.3}$$

The safe distance for the prioritized vehicle is

$$\begin{split} d\_{\text{safe,prio}}(t) &= \rho (\upsilon\_{0,\text{prio}} + a\_{0,\text{prio}}(t - \rho\_{\text{prio}})) \\ &+ \frac{(\upsilon\_{0,\text{prio}} + a\_{0,\text{prio}}(t - \rho\_{\text{prio}}))^2}{-2a\_{\text{min,break}}} - \frac{(\upsilon\_{0,\text{merge}} + a\_{0,\text{merge}}t)^2}{-2a\_{\text{max, brake}}}. \end{split} \tag{A.4}$$

The time derivatives of those terms are

$$\dot{s}\_{\text{prio}}^{+}(t) = \nu\_{0,\text{prio}} + a\_{0,\text{prio}}(t - \rho\_{\text{prio}}),\tag{A.5}$$

$$
\dot{s}^{-}\_{\text{merge}}(t) = \nu\_{0, \text{merge}} + a\_{0, \text{merge}}t \tag{A.6}
$$

and

$$\begin{split} \dot{d}\_{\text{safe,prior}}(t) &= -a\_{0,\text{prio}}\rho\_{\text{prio}} + \frac{2(\upsilon\_{0,\text{prio}} + a\_{0,\text{prio}}(t-\rho\_{\text{prio}}))a\_{0,\text{prio}}}{-2a\_{\text{min,brake}}} \\ &- \frac{2(\upsilon\_{0,\text{merge}} + a\_{0,\text{merge}}t)a\_{0,\text{merge}}}{-2a\_{\text{max,brake}}}. \end{split} \tag{A.7}$$

Setting the time derivative of the safety reserve to zero yields

$$t\_{\rm crit} = \frac{-\nu\_{\rm O,merge} + \nu\_{\rm O,prio} - \frac{\nu\_{\rm O,prio}a\_{0,prio}}{a\_{\rm min,bek}} + \frac{a\_{0,prio}^2 \rho\_{\rm prio}}{a\_{\rm min,bek}} + \frac{\nu\_{\rm O,mcye}a\_{0,mcye}}{a\_{\rm max,bek}}}{a\_{\rm O,merg} - a\_{\rm O,prio} + \frac{a\_{0,prio}^2}{a\_{\rm min,bek}} - \frac{a\_{0,mcye}^2}{a\_{\rm max,bek}}} . \tag{A.8}$$

From demanding the critical, i.e. smallest, safety reserve to be at least zero

$$d\_{\text{reserve}}(t\_{\text{crit}}) \stackrel{!}{>} 0,\tag{A.9}$$

such that no emergency maneuver is to be performed, we can compute the required initial distance of the rear vehicle *d* + init,prio, such that the merging can be performed safely when accelerating with at least *a*0,merge from the time of merging *t* = 0.

## **A.2 Evaluation Parameters**

### **A.2.1 Safety Parameters**

Table A.1: RSS Parameters.

### **A.2.2 Comfort and Convenience Parameters**

Table A.2: General Planning Parameters.

Table A.3: Accepted deceleration probabilities.

Table A.4: Parameters for the computation of vdes.

Table A.5: IDM Parameters.

Table A.6: Updated parameters for the lane change experiment, since lane changes in dense traffic require driving more dynamically.

Table A.7: Additional parameters for the narrowing experiment.

# **Bibliography**

*Intelligent Vehicles Symposium (IV)*, Dearborn, MI, USA, Jun. 2014, pp. 420–425, DOI: http://dx.doi.org/10.1109/IVS.2014. 6856487.

*Symposium (IV)*, Seoul, South Korea, Jun. 2015, pp. 322–329, DOI: http://dx.doi.org/10.1109/IVS.2015.7225706.

*Journal of Robotics Research*, vol. 26, no. 2, pp. 141–166, Feb. 2007, DOI: http://dx.doi.org/10.1177/0278364906075328.

Perception-Driven Autonomous Urban Vehicle," in *Journal of Field Robotics*, vol. 25, no. 10, pp. 727–774, Oct. 2008, DOI: http://dx.doi.org/10.1002/rob.20262.

Framework / Fusion von Distanzmessungen und semantischen Größen im Rahmen der Evidenztheorie," in *tm - Technisches Messen*, vol. 86, no. s1, pp. 102–106, Sep. 2019, DOI: http: //dx.doi.org/10.1515/teme-2019-0052.

available via: https://papers.nips.cc/paper/4031-monte-carlo-pl anning-in-large-pomdps.pdf, last retrieved 2020-05-28.

2020, *arXiv:2002.01254v1*, https://arxiv.org/abs/2002.01254v1, last retrieved 2020-05-28.

"Making Bertha Drive—An Autonomous Journey on a Historic Route," in *IEEE Intelligent Transportation Systems Magazine*, vol. 6, no. 2, pp. 8–20, 2014, DOI: http://dx.doi.org/10.1109/ MITS.2014.2306552.

# **Publications by the Author**

*Transportation Systems (ITSC)*, Maui, HI, USA, Nov. 2018, pp. 575–582, DOI: http://dx.doi.org/10.1109/ITSC.2018.8569658.

# **Supervised Theses**

Chenyang Zhou, "Decentralized Planning of Macro-Actions for Cooperative Automated Vehicles with Hierarchical Monte Carlo Tree Search," Master's Thesis, Karlsruhe Institute of Technology, Institute of Measurement and Control, Karlsruhe, Germany, Feb. 2018. Supervised jointly with Karl Kurzer.

Florian Engelhorn, "Dezentrale kooperative Planung autonomer Fahrzeuge auf Basis einer kontinuierlichen Monte Carlo Tree Search," Master's Thesis, Karlsruhe Institute of Technology, Institute of Measurement and Control, Karlsruhe, Germany, Mar. 2018. Supervised jointly with Karl Kurzer.

Bernhard Trünkle, "Entwurf und Parameteroptimierung eines Kostenfunktionals für die Trajektorienplanung von automatisierten Fahrzeugen," Master's Thesis, Karlsruhe Institute of Technology, Institute of Measurement and Control, Karlsruhe, Germany, Aug. 2018.

Rudolf Ehle, "Entwurf eines Sampling-Verfahrens für die Trajektorienplanung automatischer Fahrzeuge mittels Simulated Annealing," Master's Thesis, Karlsruhe Institute of Technology, Institute of Measurement and Control, Karlsruhe, Germany, Nov. 2018.

Gunnar Diestmann, "Trajektorienplanung für einen kooperativen Fahrstreifenwechsel hochautomatisierter Fahrzeuge," Master's Thesis, Karlsruhe Institute of Technology, Institute of Measurement and Control, Karlsruhe, Germany, May 2019.

# **Schriftenreihe Institut für Mess- und Regelungstechnik Karlsruher Institut für Technologie (1613-4214)**

Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.

Die Bände sind unter www.ksp.kit.edu als PDF frei verfügbar oder als Druckausgabe bestellbar.

Band 022 Loose, Heidi Dreidimensionale Straßenmodelle für Fahrerassistenzsysteme auf Landstraßen. 2013 ISBN 978-3-86644-942-8

### Schriftenreihe

INSTITUT FÜR MESS- UND REGELUNGSTECHNIK KARLSRUHER INSTITUT FÜR TECHNOLOGIE (KIT)

### 049

This work targets the problem of motion planning for automated vehicles. As a prerequisite for their on-road deployment, automated vehicles must show an appropriate and reliable driving behavior in mixed traffi c, i.e. alongside human drivers. Besides the uncertainties resulting from imperfect perception, occlusions and limited sensor range, also the uncertainties in the behavior of other traffi c participants have to be considered.

Related approaches often employ a deterministic problem formulation, where large uncertainties lead to conservative and over-cautious behavior. On the other hand, due to the need for real-time capability, a holistic uncertainty treatment in probabilistic approaches often induces a strong limitation of the action space of automated vehicles. Moreover, safety and traffi c rule compliance are often not considered.

Thus, in this work, three motion planning approaches and a scenariobased safety approach are presented. The safety approach enhances an existing concept - which targets the guarantee that automated vehicles will never cause accidents - by the consideration of traffi c rules for crossing and merging traffi c, occlusions, limited sensor range and lane changes. The three presented motion planning approaches are targeted towards the different predominant uncertainties in different scenarios, while operating in a continuous action space. The performance of the presented approaches is shown in various scenarios.

ISSN 1613-4214 ISBN 978-3-7315-1070-3

Gedruckt auf FSC-zertifi ziertem Papier 9<PHONE_NUMBER>03

Maximilian Naumann Probabilistic Motion Planning for Automated Vehicles


**6. Benign meningiomas or WHO grade I meningiomas**

In most of cases, WHO grade I meningiomas are diagnosed. They present a benign clinical course, with rare mitoses and occasional pleomorphic nuclei. Various architectural patterns have been described, whose most frequent are meningothelial, fibroblastic, and transitional meningiomas. Meningothelial or syncytial meningiomas are characterized by a high cellular density with cells disposed in packed sheets. Fibroblastic meningiomas are similar to schwan‐ nomas, with elongated cells. Cellular whorls and psammoma bodies (mineralized whorls) are most common in the transitional type, which present intermediate features between the meningothelial and the fibroblastic subtype. Other rare variants are the psammomatous meningiomas, characterized by multiple psammoma bodies, the angiomatous meningiomas, which are extremely vascularized with scarce foci of meningothelial cells and the secretory meningiomas, with lumina formed by cytokeratin-immunoreactive cells and containing PAS positive material.

Benign meningiomas are considered cured, if a complete resection is performed. Thus, a simple surveillance is indicated in cases where a gross total resection is possible (Simpson grade 1 or 2) with serial cerebral MRI, without any complementary treatment. A very small portion of benign meningiomas progresses to more aggressive variants (<2%) [38], but when consider‐ ing only recurrent tumors, the prevalence becomes significantly higher (14%) [39]. The clinical behavior and the malignant potential of meningioma are still matter of debate. Many au‐ thors applied the model of clonal progression to meningiomas with a progressive appear‐ ance of more aggressive cellular subpopulations according with the presence of genetic imbalances and molecular modifications [40, 41].

Meningiomas show immunohistochemical positivity for epithelial membrane antigen (EMA) in 80% of cases. They also express other epithelial markers such as vimentin and cytokeratin. Antileu 7 and glial fibrillary acidic protein (GFAP) stains, characteristics for schwannomas and gliomas, respectively, are uniformly negative. Immunohistochemistry may thus help in the differential diagnosis in difficult cases.

An elevated labeling for Ki-67 may denote a more aggressive lesion [42, 43].