Ali Hamdan
email: ali.hamdan@hds.utc.fr

Reine Talj

Véronique Cherfaoui

Transition Management Between an Autonomous Vehicle and a Real Human Driver, in a Context of Take-Over Request

Keywords: Human-machine cooperative control, transition system, cooperative driving, shared lateral control, decision-making, autonomous driving

This paper presents a new shared lateral control approach to deal with the transition management between the human driver and the autonomous system. The objective of this shared control is to realize a smooth and safe switching between the two agents steering inputs during a lane keeping maneuver. The different driving modes including the transition system are detailed. Two transition modes are defined. The human driver acts on the vehicle's lateral control. The autonomous system consists of longitudinal and lateral controller developed based on Super Twisting Sliding Mode (STSM) control approach. Then, the control authority allocation is performed using the mixing shared control that permits the fusion of two inputs via a fusion parameter. To do that, a coordinator based on a decision algorithm is developed for the driving modes and the fusion parameter determination. Finally, validation of the proposed shared control is done on the "SCANeR Studio" (OKtal) professional simulator with the human in-the-loop through the "Logitech G29" steering wheel for two different scenarios.

INTRODUCTION

According to National Highway Traffic Safety Administration (NHTSA) statistics, human leads to 90% of road accidents [START_REF] Rajamani | Vehicle Dynamics and Control[END_REF]]. In such situation, when the driver is tired or distracted, he can be assisted by an Advanced Driving Assistance System (ADAS) to prevent a critical situation. These systems reduce traffic accidents and facilitate the driver's task. This new technology of automated system starts to emerge progressively in order to automate the driving activities. There are 6 levels of automation: from 0 to 5 [START_REF] Favarò | Examining accident reports involving autonomous vehicles in california[END_REF]]. The challenge now is to realize a full autonomous vehicle (level 5) where the driver is eliminated from the control loop. However, replacing the traditional vehicles by the autonomous one needs more time to be reached. The realization of full autonomous vehicles should be evaluated step by step considering many criteria: system's robustness, road safety, ethics rules and the high cost of hardware and software, etc. On the other hand, there are many levels of autonomy where the presence of human driver is essential to navigate safely on the road (SAE level 3 and 4). In level 3, the human is the main responsible of driving tasks. He supervises the scene, monitors the vehicle's control and reacts if the automated system reaches his limits for many reasons: sensor's failures, lack of localization information, bad weather, etc. However, in the higher level of automation (level 4), the system accomplishes all driving tasks without driver's intervention. The driver can retain the control action if he is available. In addition, even if the driver is trusted during autonomous mode (auto mode), a smooth transition from auto to manual mode is needed to ensure this switching. Moreover, according to the NHTSA [START_REF] Ayoub | From manual driving to automated driving: A review of 10 years of autoui[END_REF]], the autonomous system should inform the driver about its status: (1) normal/abnormal functioning, (2) Available/not available for autonomous driving mode and (3) On/Off transition request from automated system to the driver. Thus, the transition between the manual driving and autonomous driving modes creates a new challenge for the automotive industry and researches. This transition can be from autonomous to manual driving mode and vice-versa. On the other hand, the shared control is a new aspect of control in the automotive field that ensures road's safety and enhances driving performance. During this collaborative control, the human and the autonomous controller cooperate together to accomplish a common goal simultaneously. In this work, we will treat the shared lateral control in order to realize a safe and smooth transition when switching between both agents. This transition is defined as shared control where the driver and the autonomous system are in charge of driving tasks. Many advanced studies have been developed in the literature to investigate the shared control between driver and an automated system. In [START_REF] Rath | Robust lane keeping control in automated vehicles: A driver-inthe loop approach[END_REF], a shared control approach was developed based on nonlinear vehicle-road-driver model for the lane keeping purpose. A sharing parameter was integrated in the shared design to consider the conflict between the two agents. There are two types of shared control: the blending shared control and the haptic shared control. The authors in [START_REF] Johns | Exploring shared control in automated driving[END_REF] proposed a haptic-feedback shared control between the driver and the automated system using the torque on the steering wheel. The study shows that the haptic feedback torque on steering wheel is insufficient to predict the driver's intention in the future. A two-phase haptic interface based on the human-machine interaction is developed in [START_REF] Lv | Human-machine collaboration for automated driving using an intelligent two-phase haptic interface[END_REF] to manage the authority shift from autonomous to manual mode. This haptic torque on the steering wheel aims to guide and assist the human while he is engaged in the control loop, depending on his ability. However, the shared control can cause a conflict in dangerous situations. A haptic guidance torque calculated with respect to the driver's behavior is presented in [START_REF] Boink | Understanding and reducing conflicts between driver and haptic shared control[END_REF] to mitigate the conflict between both agents. Other studies represent the shared control as transition control phase affected by the human's factors. In [START_REF] Kim | Transitions from autopilot to manual control in highly automated driving: Cognitive simulations[END_REF], the authors presented a driver's cognitive model to analyze the influence of driver's experience and workload on the transition between the two agents. The authors in de [START_REF] De Winter | Preparing drivers for dangerous situations: A critical reflection on continuous shared control[END_REF] have explored the topic of shared control in terms of its advantages and drawbacks on the driving tasks. They presented the risks that can occur during this shared control, the conflict for example, and the influence of this control strategy on the driver's skills. A planned control transition was applied in [START_REF] Holländer | Preparing drivers for planned control transitions in automated cars[END_REF] to switch from automated system to the driver by using the auditory and visual information systems. These systems help and prepare the human who is engaged in non driving related activities to be ready for regaining control. Moreover, an overview for the last ten years was done in [START_REF] Ayoub | From manual driving to automated driving: A review of 10 years of autoui[END_REF], and summarizes the different works related to the shifting from manual to autonomous driving. All these interesting research have motivated us to develop a shared lateral control to address the transition management between the human and the automated system. The shared lateral control is done for the lane keeping purpose using steer-bywire system. However, this method can be adapted to consider the tracking of a local trajectory with more complex maneuvers (overtaking, collision avoidance...). Thus, in our present work, a new shared lateral control is developed to manage the transition control in order to ensure a smooth switching and enhance road safety. Our contribution is illustrated in the development of a transition system. Two modes are defined in the transition system, that ensure the shifting from auto to manual mode and vice-versa, according to many criteria: the Driver's availability (DA), the conflict on the driver's behaviors (δ con f lict ), and the take over request (TOR). Then, a coordinator based on a decision algorithm is developed for the determination of the driving mode. Finally, two scenarios are considered to test and validate the proposed shared lateral control on "SCANeR Studio" (OKtal) professional simulator. Validation is done to show the effectiveness of this method in terms of driving performance enhancing and vehicle's stability keeping. The rest of the paper is as follows: in Section 2, the different driving modes including the transition system are given. Then, the proposed shared lateral control is detailed in Section 3. A coordinator based on a decision algorithm is introduced, then the shared lateral control authority is detailed to blend the two inputs. In Section 4, system's validation is done by describing the simulations tools and presenting the simulations results of two detailed scenarios. The validation of the proposed shared control is done with a complete nonlinear model on "SCANeR Studio" (OKtal) simulator interacted with human through the "Logitech G29" steering wheel. Finally, the conclusions and the perspectives for future work are given in Section 5.

DRIVING MODES

Human Driver

A real human driver acts on the vehicle's lateral control through a steering wheel (human in-the-loop). His steering input will be noted as δ sw-h in this paper.

Autonomous System

A robust autonomous system based on the super-twisting sliding mode (STSM) algorithm is developed to realize a lane keeping maneuver at the desired velocity. The autonomous system consists of longitudinal and lateral robust controllers. The autonomous controller is validated separately on Matlab/Simulink using the robotic formalism of nonlinear vehicle model validated on "SCANeR Studio"simulator. Let us introduce briefly the STSM theory. The main idea is to define a sliding surface, representing the desired behavior of the system, where the dynamic states are forced to reach this surface during a finite time and remain on it. Consider the second order system given as: Ẋ = f (X,t) + g(X,t)u(t)

(1) where X = [x, ẋ] T ∈ ℜ 2 is the state vector, u is the control input, and f , g are continuous functions. X * is the desired state of X with X * = [x * , ẋ * ] T ∈ ℜ 2 . The error vector is given by E = X -X * = [e, ė] T ∈ ℜ 2 where e = xx * and ė = ẋẋ * . Therefore, a sliding variable s with relative degree r = 1 w.r.t the control input, is defined as:

s = ė + k e.
(2) The second order derivative of s can be expressed as:

s(s,t) = Φ(s,t) + ξ (s,t) u(t)
(3) where Φ(s,t) and ξ (s,t) are unknown bounded functions. The goal of the Super-Twisting algorithm is to force the sliding variable s to converge to zero (s = 0) in finite time. Assume that there exist positive constants S 0 , b min , b max , C 0 , U max verifying for all X ∈ ℜ n and |s(X,t)| < S 0 :

|u(t)| ≤ U max |Φ(s,t)| < C 0 0 < b min ≤ |ξ (s,t)| ≤ b max (4)
Thus, the control input based on the Super-Twisting Sliding Mode algorithm [START_REF] Rivera | Super-twisting sliding mode in motion control systems[END_REF]] is given as:

u(t) = u 1 + u 2 u 1 = -α 1 |s| τ sign(s), τ ∈]0, 0.5] u2 = -α 2 sign(s) (5)
α 1 and α 2 are positive gains. The following conditions guarantee the finite time convergence:

   α 1 ≥ 4C 0 (b max α 2 +C 0 ) b 2 min (b min α 2 -C 0 ) α 2 > C 0 b min (6)
The convergence analysis is shown in [START_REF] Utkin | On convergence time and disturbance rejection of super-twisting control[END_REF].

The controller synthesis is based on a robotic formalism model presented in [START_REF] Chebly | Coupled longitudinal/lateral controllers for autonomous vehicles navigation, with experimental validation[END_REF], Chebly (2017) (more details in Appendix A), that represents the coupling between the longitudinal and lateral dynamics. Based on this model, we choose the two sliding variables for the longitudinal and lateral controllers as follows:

s 1 = e V x + λ x e V x , λ x > 0 s 2 = ėy + λ y e y , λ y > 0 (7)
where λ x and λ y are positive constants, and, e V x (e V x = V x -V * x ) and e y (e y = y-y * ) are the vehicle longitudinal speed error and the lateral error respectively. The sliding variables s 1 and s 2 have a relative degree equal to one w.r.t the inputs respectively, the driving/braking torque Γ c for the longitudinal dynamics and the steering wheel angle δ sw-as for the lateral dynamics. Thus, in order to converge these variables to zero and the controlled states follow the desired ones, and based on the above discussion, the torque and the steering angle control applied to the vehicle, are given by:

Γ c = -α Γ c ,1 |s 1 | τ Γc sign(s 1 ) -α Γ c ,2 t 0 sign(s 1 )dτ, δ sw-as = u 1 + u 2 + δ *    u 1 = -α δ ,1 |s 2 | τ δ sign(s 2 ), u 2 = -α δ ,2 t 0 sign(s 2 )dτ, δ * is the equivalent control input, (8)
where α δ ,i and α Γ c ,i with i = [1, 2], are positive constants satisfying the conditions in (6). τ Γ c and τ δ are constants in ]0, 0.5]. The controller gains are given in [START_REF] Chebly | Coupled longitudinal/lateral controllers for autonomous vehicles navigation, with experimental validation[END_REF]. Finally, δ * is the equivalent control input, corresponding to the steering wheels angle at the equilibrium when ṡ2 = 0.

Transition System

The transition system is defined in the way to switch from the auto to manual mode or from manual to auto mode. For that, two modes are defined to realize these transitions. The transition structure is given in the Fig. 2(a). The two modes are given as:

• Transition mode 1: This mode is defined to switch from the auto mode to the manual mode after a take over request demanded by the driver (Fig. 2(b)). We define α des function of the driver's availability (DA) and the conflict on the driver's behaviors (δ con f lict ). The role of α des is to determine the direction of transition (from auto to manual mode or vice-versa, auto mode corresponds to 0, and manual mode corresponds to 1), depending on the driver's availability and behaviors. α des is the desired reference for the fusion parameter α, given as:

α des = DA * (1 -δ con f lict ) (9)
This transition is done, if the driver demand a take over request (T OR = 1) and α des = 1, that means the driver is available to take the action of vehicle's control (DA = 1), and there is no conflict on his behavior (δ con f lict = 0). So, the transition mode 1 is activated in order to switch to the manual mode and α increases from 0 to 1 in T up = 1.5s. α is the fusion parameter varying between 0 and 1 to blend the two agents steering inputs. α is equal to 0 in auto mode and to 1 in manual mode. The driver's availability (DA) is a dynamic variable related to the driver, and can be calculated based on different factors: driver's eyes movement, driver's head position, level of driver's sleepiness, etc. Therefore, the calculation of diver's availability is not in the scope of this work and it is considered as an input to the transition system. On the other hand, δ con f lict is the difference between the steering wheel angle applied by the driver (δ sw-h ) and the approximate steering angle (δ approx ) on the desired trajectory given as:

δ con f lict = B δ threshold (|δ approx -δ sw-h |) (10)
δ approx = ρ * * (l f + l r ) (11)
with ρ * the curvature of the desired trajectory, l f and l r the distances from the center of the vehicle to the front and rear wheels respectively. The conflict of the driver's behavior, δ con f lict is a Boolean value, equal to 0 or 1 depending on δ threshold a threshold from which we consider that there is a conflict on the driver's behavior. • Transition mode 2: The aim of this mode is to assure the transition from the manual mode to the auto mode, in the case when the human driver is totally engaged in the driving tasks (manual mode) and suddenly the value of α des is equal to 0 (Fig. 2(b)). So, the transition mode 2 is activated to realize the transition from manual mode to auto mode and α decreases from 1 to 0 in T down = 0.2s. In this case, the transition mode 2 is an example of Advanced driving Assistance system (ADAS) application.

In addition, the transition mode 2 is activated if during the transition from auto to manual mode (transition mode 1), suddenly the value of α des is equal to 0, that means if a conflict occurs on the driver's behaviors, or the driver is not available during this phase (Fig. 2(b)). Thus, the transition to the auto mode is done and α decreases from α 0 (α 0 < 1) to 0, respecting the same decrease speed.

SHARED LATERAL CONTROL ARCHITECTURE

This section details the proposed shared lateral control (see Fig. 1). A coordinator based on a decision algorithm is developed to determine the way of switching between the different driving modes. Then the control authority between the two agents is presented.

Coordinator

In order to determine the driving mode and the value of the fusion parameter, a coordinator is needed to switch between the different driving modes: Auto mode, manual mode and finally the modes 1 and 2 of the transition system (see Fig. 3). These modes are defined as:

Auto mode

Manual mode Mode 1

Mode 2

Coordinator

Transition system Driving modes α Fig. 3. The structure of the coordinator • Auto mode: an autonomous controller provides the vehicle with the appropriate torque for the longitudinal movement. For the lateral dynamics, a robust controller generates the steering angles to follow the desired trajectory and keep the lane. Noting that, in this mode the steering wheel and the wheels are decoupled corresponding to a steer-by-wire system. • Manual mode: The human driver acts on the vehicle's lateral control by using the steering wheel, "Logitech G29" (Fig. 6(a)), while the longitudinal controller is realized by the autonomous controller. The driver can demand a take over action by pressing the bottom R2 and he can turn off the manual driving by using the bottom L2 (see Fig. 6(b)).

• Transition system's modes: transition mode 1 and 2 are detailed above to describe the switch from Auto mode to the Manual mode and vice-versa.

Then, a decision algorithm (given below) is developed in this layer for the switching between the different modes, and the determination of the final value of fusion parameter α.

Shared Lateral Control Authority

The shared control authority allocation between the driver and the autonomous system is performed by using the blended shared control [START_REF] Borroni | A weighting approach to the shared-control of lateral vehicle dynamics[END_REF]] which allows the fusion of two inputs from each agent. This form of shared control is used in the steer-by-wire system, where there is no mechanical link between the two inputs. The familiar form of blending shared control is the mixing using weight parameter (see. [START_REF] Li | Shared steering control of tire blowout for ground vehicles[END_REF]], [START_REF] Li | Shared control driver assistance system based on driving intention and situation assessment[END_REF]). The total blending control input is given as:

δ total = α * δ sw-h + (1 -α) * δ sw-as , (12
) where δ total , δ sw-h and δ sw-as are the total steering wheel angle, steering wheel angles of human driver and autonomous system respectively. α is the fusion parameter representing the Algorithm 1 A decision algorithm for the switching purpose if (Take over request equal to 0 (TOR=0)) then Auto mode is activated

(α = 0) else if (α des = 0) then if (α t-1 = 0) then Auto mode is activated else Transition mode 2 is activated end if else if (α t-1 < 1) then Transition mode 1 is activated else Manual mode is activated (α = 1) end if end if end if
influence proportion of each agent on the total steering angle. α is bounded in [0,1], calculated depending on the decision algorithm in the coordinator.

SYSTEM'S VALIDATION

The simulation tools used in the validation are presented in this section. Validation is done on the"SCANeR Studio" (OKtal) simulator by co-simulation between the later and Matlab/Simulink. Then, the simulation results of two scenarios are presented later to show the effectiveness of the proposed shared lateral control including the transition system.

Simulation tools

-SCANeR Studio simulator: ''SCANeR Studio" simulator is a simulation platform answering the demand of researches and engineers. This simulator allows user to create a safe and controllable environment for the validation of different scenarios with the different driving conditions (Fig. 4). A Scaner vehicle model is used on this software, with the different driving modes: Autonomous and Manual modes and the 2 modes of transition system. Thus, The "SCANeR Studio" simulator is used in this work to validate the proposed shared control that interacts with the human driver through the "Logitech G29" (see Fig. 6).

-Matlab/Simulink:

The co-simulation is done on Matlab/Simulink, by developing the different Simulink blocks in Fig. 5. The "Logitech G29" (Fig. 6) is the steering wheel that permits the interaction between driver input with the "SCANeR Studio" simulator. The wheel can rotate up to 900 • , corresponding to 60 • on the wheels. It has three pedals and a dual motor force feedback with a overheat safe guard. All these features make the driving more realistic for the users. In addition, the different blocks for the autonomous controller, localization, etc., are developed on Matlab/Simulink.

Simulation results

After preparing the simulation environment to validate the proposed shared control, validation is done on "SCANeR Studio" simulator, by using the test track given in the Fig. 7. The vehicle realizes the path following maneuver by cooperating with the human in two different scenarios discussed in the following. Note that δ threshold =1.2 rad for the 2 scenarios. -Scenario 1: This scenario shows the transition between the auto and manual modes depending on the driving situations. The transition is able to switch from one to another mode by ensuring road safety and enhancing driving performance. The driver's availability is full for the overall trajectory (see Fig. 8). The driver demands a take over request (TOR=1 at t=8.5s) by pressing the R2 bottom and α des = 1 because the driver's availability is equal to 1 and there is no conflict on the driver's behaviors (eq.9). For that, α increases from 0 to 1 in 1.5s at t=10s according to the transition mode 1 of transition system. The driver takes the action of vehicle's control until t=32s. At t=32s, there is a conflict on driver's behaviors (see. Fig. 9) and α des becomes 0. For that, α decreases again to 0 in 0.2s to penalize the driver, and the driving will be fully autonomous, according to the transition mode 2. Again at t=50s the driver still asking a take over action and α = 1 because the conditions of activation of transition mode 1 are realized. The driver acts on the vehicle's control until t=70s and finally he decides to switch off the manual mode by pressing the L2 bottom and α=0 at t=70s. The different values of α are given in the Fig. 8, which shows that the lateral error is between -40cm and +40cm, that means the cooperation of driver with the autonomous system leads to an acceptable and accurate lane keeping with a small lateral error. The steering angles on the wheels of Auto, Manual and shared modes are given in the Fig. 9. As we can see in Fig. 9, a smooth switching is done between the two agents to ensure the transition between them and compensate the conflict of driver's behaviors, that makes the system stable. Finally, the Fig. 10 shows the longitudinal speed which tracks the desired one using the STSM controller, the road curvature of the desired track and finally the lateral and longitudinal accelerations. The actual lateral acceleration does not exceed the ±5m/s 2 , which corresponds to a comfortable driving zone.

In addition, the actual longitudinal acceleration is pertinent (< ±3m/s 2 ) for a comfortable maneuver. To conclude, the Scenario 1 proves that the transition between the two agents is possible while maintaining a comfortable driving against the switching between the different modes and enhancing safety on the road.

-Scenario 2: This scenario is defined to deal with the variation of different variables (DA, δ con f lict , α des , TOR). The transition system is adapted against this variation to show the functionality of this system. The management between the different modes under varying driving conditions is presented here. Refer to Fig. 11, the driver asks the take over action permission at t=5s, however this action is not approved until t=10s where the driver's availability is full (see Fig. 11) and α des = 1. α starts to increase and the driver takes totally control at t=11.5s. He is still acting on the vehicle's control until t=25s, where the DA is lost (α des = 0) and the auto mode is switched on to drive the vehicle until t=40s. At t=40s, the driver retains the vehicle's control because the DA is equal to 1 again (α des = 1) and α starts to increase. During this phase, there is a conflict that occurs on the driver's behaviors (α des becomes 0) and the value 1 for α is not reached. α will decrease again to 0 to give the control authority to the autonomous system in case the driver is not able to drive the vehicle. For the same reasons, α has the same behavior around t=50s. In addition, α=0 at t=60s even though the conflict on the driver's behaviors starts to decrease (see Fig. 12), but the driver's availability is equal to 0 again and α des = 0 (eq.9) at t=60s. The increase of DA at t=70s leads to increasing in the α value. The driver drives his vehicle from t=70s to t=83s when he decided to switch on the Auto mode (TOR=0) and take his rest. The parameter α is given in the Fig. 11. As shown in Fig. 11, the different transitions lead to a small lateral error that guarantees the driving safety in terms of vehicle's stability during the lane keeping. Fig. 12 shows the different steering angles on the wheels. The transition between the manual and auto modes is smooth and safe specially during the switching (at t=25s, 40s, 70s for example). As we can see in the Fig. 13, there is a good tracking for the desired longitudinal speed. The actual lateral and longitudinal acceleration are limited between ±5m/s 2 and ±3m/s 2 respectively, which confirm the passenger's comfort. Thus, the driving is safe and comfortable despite the variation of driving conditions and the switching between the three modes: auto, manual and transition system's modes.

CONCLUSION AND PERSPECTIVES

To conclude, in this paper a shared lateral control has been developed to ensure the transition between the human driver and the automated system. The different driving modes including the transition system are detailed in this work. Then, the shared lateral control architecture is presented to blend the two control inputs. A decision making algorithm is developed, using a coordinator for the control authority allocation and driving modes determination. The proposed shared control is validated on the "SCANeR Studio" (OKtal) professional simulator interacting with human driver through the "Logitech G29" steering wheel. The validation is done for the two scenario of lane keeping maneuver. The results show the effectiveness of the proposed approach to ensure a smooth transition and promote road's safety. In the future work, we will consider other criteria to prove the effectiveness of the decision making process during the transition phase.

Fig. 1 Fig. 1 .

 11 Fig.1shows the global shared lateral control architecture. It is composed of two parts: the Driving modes and the shared lateral control. The different driving modes are detailed in the following.

Fig. 2 .

 2 Fig. 2. The structure of the transition system

Fig. 4 .

 4 Fig. 4. The SCANeR Studio environment.

Fig. 6 .

 6 Fig. 6. The steering wheel hardware, Logitech G29.

Fig. 7 .

 7 Fig. 7. Map of the test track.

Fig. 13 .

 13 Fig. 8. The lateral error -S1

ψ) 2 )δ . The control inputs to the vehicle are the steering wheel angle, δ , and the Driving/Braking wheels torque τ w . More details about the vehicle parameters are given in[START_REF] Chebly | Coupled longitudinal/lateral controllers for autonomous vehicles navigation, with experimental validation[END_REF].

Appendix A. REDUCED VEHICLE MODEL

The vehicle model is based on a robotic formalism model presented in [START_REF] Chebly | Coupled longitudinal/lateral controllers for autonomous vehicles navigation, with experimental validation[END_REF], [START_REF] Chebly | Trajectory planning and tracking for autonomous vehicles navigation[END_REF]. However, to simplify the controller design, some assumptions are given in [START_REF] Chebly | Coupled longitudinal/lateral controllers for autonomous vehicles navigation, with experimental validation[END_REF]. With all assumptions, the reduced vehicle model used to establish control laws, can be written as:

where x and y are the longitudinal and the lateral positions of the vehicle at its center of gravity (COG). ψ is the yaw angle.

F aero is the longitudinal aerodynamic force. m e , g 1 , g 2 and g 3 are given by:


Bart M H Bruininks
email: b.m.h.bruininks@rug.nl

Albert S Thie

Paulo C T Souza

Tsjerk A Wassenaar

Shirin Faraji

Siewert J Marrink

Sequential voxel-based leaflet segmentation of complex lipid morphologies

As molecular dynamics simulations increase in complexity, new analysis tools are necessary to facilitate interpreting the results. Lipids, for instance, are known to form many complicated morphologies due to their amphipathic nature, becoming more intricate as particle count goes up. A few lipids might form a micelle, where aggregation of tens of thousands could lead to vesicle formation. Millions of lipids make up a cell and its organelle membranes, and are involved in processes such as neurotransmission and transfection. To study such phenomena, it is useful to have analysis tools which understand what is meant by emerging entities such as micelles and vesicles. Studying such systems at the particle level only, becomes extremely tedious, counter-intuitive and computationally expensive. To address this issue we developed a method to track all the individual lipid leaflets, allowing for easy and quick detection of topological 1 changes at the meso scale. By using a voxel-based approach and focussing on locality, we forgo costly geometrical operations without losing important details and chronologically identify the lipid segments using the Jaccard index. Thus we achieve a consistent sequential segmentation on a wide variety of (lipid) systems, including monolayers, bilayers, vesicles, inverted hexagonal phases, up to the membranes of a full mitochondrion. It also discriminates between adhesion and fusion of leaflets. We show that our method produces consistent results without the need for pre fitting parameters, and segmentation of millions of particles can be achieved on a desktop machine.

Introduction

The increase in complexity of molecular dynamics (MD) simulations poses challenges for current analysis software. This is well illustrated in the field of biological applications. Such simulations changed from a single protein in vacuum 1 in the late 1970s to complex membrane based systems containing millions of molecules such as (phospho)lipids, cholesterol, proteins, DNA, drugs, water, etc. [START_REF] Marrink | Computational Modeling of Realistic Cell Membranes[END_REF][START_REF] Van Eerden | Characterization of thylakoid lipid membranes from cyanobacteria and higher plants by molecular dynamics simulations[END_REF][START_REF] Bruininks | A molecular view on the escape of lipoplexed DNA from the endosome[END_REF][START_REF] Vögele | Membrane perforation by the pore-forming toxin pneumolysin[END_REF][START_REF] Ingolfsson | Lipid Organization of the Plasma Membrane[END_REF][START_REF] Perilla | All-Atom Molecular Dynamics of Virus Capsids as Drug Targets[END_REF] The number of molecular aggregates in simulations has also increased, and nowadays it is not rare to see a mixture of micelles, bilayers and vesicles in a single simulation. [START_REF] Bruininks | A molecular view on the escape of lipoplexed DNA from the endosome[END_REF][START_REF] Ramezanpour | Structural Properties of Inverted Hexagonal Phase: A Hybrid Computational and Experimental Approach[END_REF][START_REF] Holdbrook | Dynamics of Crowded Vesicles: Local and Global Responses to Membrane Composition[END_REF][START_REF] Pezeshkian | Backmapping triangulated surfaces to coarse-grained membrane models[END_REF][START_REF] Alasiri | Determining Critical Micelle Concentrations of Surfactants Based on Viscosity Calculations from Coarse-Grained Molecular Dynamics Simulations[END_REF] This increase in complexity of both composition and aggregation state requires a higher level of abstraction in analysis than typically offered. During analysis we would like to work with concepts such as bilayers, vesicles, leaflets, protein aggregates, phases, etc. on top of simple identifiers such as molecule names and indices, often used in current tools. This would dramatically reduce the complexity of coding the analysis software and would allow for a more user friendly interface with data. Current MD analysis software, such as MDAnalysis, [START_REF] Gowers | A Python Package for the Rapid Analysis of Molecular Dynamics Simulations[END_REF] Pytim, [START_REF] Sega | Pytim: A python package for the interfacial analysis of molecular simulations[END_REF] FATSLiM 14 and VMD, [START_REF] Humphrey | Visual molecular dynamics[END_REF] come with built-in functions or have packages to detect e.g. bilayer leaflets and aggregates, and offer some of the desired level of abstraction. However, these tools were developed to deal with relatively simple compositions and topologies, and may fail as systems contain more than one or two segments (SI Fig. 8).

To avoid this problem we developed a robust method for consistent sequential segmentation with a voxel oriented analysis paradigm. No assumptions are required regarding the number of segments or molecules, and it natively supports all-atom or coarse-grained (CG) simulation data. Furthermore, voxel-based methods offer very good scaling with the increase of particles and are excellent for locality queries. [START_REF] Dado | Geometry and Attribute Compression for Voxel Scenes[END_REF][START_REF] Villanueva | Proceedings of the 20th ACM SIGGRAPH Symposium on Interactive 3D Graphics and Games[END_REF][START_REF] Careil | Interactively Modifying Compressed Sparse Voxel Representations[END_REF][START_REF] Sualeh | Dynamic Multi-LiDAR Based Multiple Object Detection and Tracking[END_REF] To guide the reader through this work we will first present pseudo code to illustrate the general implementation of our leaflet segmentation algorithm, which is implemented in a tool coined MDVoxelSegmentation. We then show a wide range of examples, and we end with a discussion of the limitations and further prospects of the tool. The Supporting Information (SI) contains the video results of the algorithm.

Design and implementation

In the following sections, we describe the complete algorithm in two sections of pseudo code.

We start with considering a single frame and perform our multi-step connective components analysis for lipid-based systems. Then we describe how we add the extra layer of chronological consistency on top.

Pseudo code for spatial segmentation of leaflets

The aim is to segment a lipid density into leaflets (i.e. a single membrane leaflet is a segment).

To successfully segment lipid densities into leaflets we require three things, a reference, a trajectory and a selection file. The reference file should contain information regarding atoms and their naming (e.g. GRO, PDB). The trajectory file should contain coordinates for each point in the reference file in every frame (e.g. GRO, XTC, TRR). Finally, the selection file should contain the selections (lipid heads, lipid tails, lipid linkers and exclusion beads). By default most CG Martini lipids are supported, [START_REF] Wassenaar | Computational Lipidomics with insane: A Versatile Tool for Generating Custom Membranes for Molecular Simulations[END_REF] meaning no selection has to be specified for the leaflet segmentation of those systems. However, if you are working on a custom lipid, it can easily be added to the default selections file (an example for CHARMM lipids is given in the methods for the leaflet segmentation). The completed segmentation will assign each lipid (residue) a segment identity. A schematic representation of leaflet segmentation is shown in Fig. 1.

Leaflet segmentation starts with reading the input files (Fig. 1 A) and mapping all the lipid tails, lipid heads, and exclusion beads to boolean voxel matrices (0.5 nm binning for Martini lipid systems; Fig. 1 B). This results in the following three 3D boolean matrices: tails, heads and exclusions. The exclusions are always grown once, meaning all neighboring voxels (26) of a True voxel become True. The exclusions are used as local stops during segmentation, acting as boundaries for the segmentation (i.e. exclusions can be used to confine the segmentation). We found that using exclusions can be beneficial when working with proteins, since the protein boundaries are prone to be edge cases for lipid segmentation and are better handled later.

The tails matrix is segmented using the adapted connected components algorithm (Fig. 1 C). Before segmentation starts, the heads matrix is subtracted from the tails matrix. The exclusions act as boundaries during segmentation. The result of the tails segmentation is depicted in Fig. 1 D.

Roughly the same procedure is followed for segmenting the heads as was done for the tails. However, instead of using the complete heads matrix, it is first deconstructed by the tails segments. Each tails segment is individually processed (Fig. 1 E,G). This is achieved by using only those headgroups which have a lipid residue index which is also in the active tails segment. Before heads segmentation starts, the tails matrix is subtracted and the exclusions act as boundaries during segmentation. Segmenting the heads in this manner results in the leaflet segmentation per tails segment (Fig. 1 F,H).

Although the result at this point is good (Fig. 1 I,J), there might be some lipids which still have no segment assigned. Force-segmentation can solve this problem as it can be used to map lipids to a segment based on a distance criterion. For every lipid which does not have a segment assigned, all neighbouring lipids are checked within a certain cutoff distance (default is 2 nm). If the dominant segment in that search is not 0 (unsegmented), the lipid is appointed to that dominant segment. Our implementation performs this operation in an iterative manner by trying to classify with small cutoff distances first (default is 1 nm).

The cutoff is slightly increased (default is 0.1 nm) if no extra lipids were segmented. If an iteration does result in new lipids being segmented the cutoff is set back to the minimum value. The iteration halts if the force-segmentation returns a stable output (i.e. no further changes occur), if the cutoff becomes larger than the maximum value, or if all selected lipids have been segmented. The final result after force-segmentation is depicted in Fig. (1 K).

Force-segmentation can be turned off.

Splitting the segmentation into two steps, using force-segmentation as a post processing step, means that simple cases are handled in a relatively low-cost manner where hard cases might use a different approach altogether (in this case our force-segmentation protocol). This mixed use of a low-cost algorithm with more precise and expensive algorithms reduces the total complexity of the calculation, whilst at the same time hardly affecting the segmentation quality.

The result after spatial leaflet segmentation of a trajectory, is a new trajectory which contains a segmentation entry for each index in each frame. These indexes correspond to the atom indices in the original structure file. Atoms which have no segmentation are assigned to 0. However, due to the arbitrary picking of our initial voxels for segmentation, the segmentation labels miss sensible consistency over time for all other labels than 0. To solve this lack of chronological consistency, we need to perform temporal segmentation.

Pseudo code for temporal segmentation

All that is required as input for the temporal segmentation is the spatial segmentation output.

To perform temporal segmentation, we need to find the identity of segments for each frame sequentially. When comparing a new frame with a current frame, the algorithm first considers the identities of segments of the current frame. The identified segments of the current frame are then matched to the segments in the new frame (Fig. 2 A). The atoms within a segment are used as the basis for the identity of the segment. Segments are treated as sets, for which the unique elements are the atom indices. The identity of a segment in the new frame is the segment in the current frame with which it shares the most atoms. In this way segments can be matched across frames, even if due to shrinking or growth their atom contents are not identical. Consider frame N which contains two segments N 1 = {a1, a2, a3, a4} and N 2 = {a5, a6, a7, a8} and frame M which contains the segments M 1 = {a1, a2, a3, a8} and M 2 = {a5, a6, a7, a4, a9}. Each segment contains unique elements, which have a lowercase unique code. When we compare frame N with frame M, we see that the elements a8 and a4 have swapped places and a9 is introduced to M2. However we still want to assign N1 to M1 and N2 to M2, to allow for the growing and shrinking of segments. We therefore need a way to identify the amount of shared elements between two segments.

J(A, B) = A ∩ B A ∪ B (1)
To identify how many atoms are shared between two segments across frames the Jaccard index is used (eq. 1). This index is defined as the intersection of two sets divided by the union of two sets. This provides a ratio between 0 and 1, where 1 is complete similarity and 0 means no shared atoms. The parameter Jaccard threshold defines the minimum score on the Jaccard index that is required for two sets to have the same identity (SI Fig. 10). For example a Jaccard threshold of 0.618 means that only if two segments have a score higher than 0.618 together, they have the same identity. This allows for consistency while also identifying large shifts between frames in the simulation. For example a leaflet in a bilayer system which has lost lipids due to flip-flopping can still be indentified in the next frame.

Because of the possibility of these large shifts the algorithm does not assume that the amount of segments between frames remains the same. It therefore allows for the disappearance of identified segments and the appearance of new identities. Identified segments can disappear when no match in the new frame can be made that exceeds the Jaccard threshold. An identified segment disappearing is considered to have merged with another segment. A new identity can be created when in the new frame a segment exists to which no old identified segment has matched. This new segment then receives a new identity. A newly created identified segment is considered to have split from a previously existing identified segment. In the case of a merge (Fig. 2 B), two or more identified segments in the current frame are merged into an unidentified segment in the new frame. Two distinct cases are possible in the new frame. Either one of the identities of the merging segments is inherited by the new segment or a new identity is created. If one of the identities is maintained, the Jaccard score between one of the identified segments and the unidentified segment has exceeded the Jaccard threshold. This indicates a large segment is merging with one or more smaller segments. For example, this could be the result of a transient fusion between two leaflets.

If a new identity is created, this means none of the previous segments contained enough atoms to attain a Jaccard score higher than the Jaccard threshold when compared to the unidentified segment. In this case all previous identified segments are merged into a new segment with a new identity.

To allow for the reappearance of identities, the disappeared identities are stored in a database. This database contains the identity and the atoms present in the identified segment at the current frame. It also logs into which segment the disappeared segment has merged.

The segment into which the disappeared segment has merged is the segment for which it had the highest Jaccard score. This score will always be below the Jaccard threshold, as otherwise the identity would be inherited. This information can be used to restore identities when splitting events take place.

In a split, two or more unidentified segments in the new frame have the highest Jaccard score with the same identified segment in the current frame (Fig. 2 C). In such a split there are also two distinct cases. In one case a small part splits off from an existing identified segment. However one of the unidentified segments still has a high enough Jaccard score to match with the identified segment. This means that one of the unidentified segments inherits the identity of the identified segment. This could indicate a leaflet loosing a small group of lipids which form their own segment. The other unidentified segments will then either gain a new identity or take the identity of an old segment from the database. The new segment takes the identity of the segment with the highest Jaccard score in the database, but only if the score exceeds the Jaccard threshold. If the threshold is not exceeded, the segment is assigned a new identity.

In the other case the identified segment could not match with one of the unidentified segments. This will lead to the disappearance of the identified segment. This could be caused by a leaflet splitting into two roughly equal parts. For example, in the case of a bilayer with a transient pore, closing of the pore would result in a split of the bilayer in its previous two leaflet segments. The disappearing segment will be stored in the database, while the new segments will each be compared to the database. Each segment which has a match with the database will continue with the identity of this match, while the unmatched segments continue with new unique identities.

The result after temporal segmentation is a consistent sequential segmentation over the trajectory, both spatially and over time. Identities between time frames are consistent, based on the atoms within each segment. Merging segments are detected and identified in a size dependent way. Splitting segments are detected as well. The database of previous segments allows identities to persist over a longer time frame.

Results

To validate the leaflet segmentation algorithm, we prepared a wide range of lipid systems using the Martini CG force field. [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF] The same default input was used for all segmentation unless specifically stated otherwise. We start with simple lipid bilayer systems in varying topologies (Fig. 3). Then we add proteins and a high variety of compositional complexity (Fig. 4). This is followed by a complex phase transition and fusion of DNA-lipid complexes (lipoplexes) with endosome membrane models to illustrate segmentation quality in the presence of tight, curved and complex dynamic geometry (Fig. 5). We end with some more general real world examples, including an example at atomistic resolution (Fig. 6).

Pure bilayer systems

Leaflet assignment for a simple bilayer containing 270 DOPC and 66 cholesterol (4:1) without force-segmentation resulted in good segmentation, but some lipids, in particular cholesterol, are unlabeled (Fig. 3 A, C orange). By turning on force-segmentation all lipids including the cholesterol could be assigned successfully to a leaflet (Fig. 3 B, C green).

A common lipid property to investigate is the flip-flop rate. We compared the assigned flip-flopping of our leaflet segmentation to a common procedure. This common analysis would assign a lipid flip-flop based on the z-coordinate of a lipid headgroup with respect to the average z-position of the bilayer. To prevent high frequency noise due to headgroups at or close to the membrane center, a deadzone is used. Lipids within the deadzone region around the center are not assigned to a leaflet. The downside of the z-height approach is that it can only be used for small patches of a relatively flat bilayer with its normal pointing along the z dimension. Our method takes a rather different approach. We check for changes in leaflet membership for each lipid. To reduce noise, a lipid which is not segmented is assigned to its previous valid leaflet. The advantage of our segmentation based method over the conventional approach is that it still works if the membranes/leaflets are curved, or if they are not a well defined bilayer at all. In our comparison of cholesterol flip-flopping, our algorithm performs well which is indicated by the close relation between the segmentation labels and the discretization of the lipid z-height (Fig. 3 C). Segmentation without forcesegmentation is good enough for lipid flip-flop calculations. The segmentation based flip-flop analysis compares well with the z-height approach if a deadzone of ±0.6 nm is used in the latter (Fig. 3 D). Although the flip-flop patterns for both methods are highly similar, they are not the same, therefore, direct quantitative comparison between both methods should be performed with care. However, if a z-height deadzone of ±0.3 nm is used, nearly all flip-flops indicated by the segmentation analysis are in the set of the z-height approach (SI Fig. 9).

In other words it appears that the flip-flops indicated by the segmentation are a near perfect subset of the z-height approach with a relatively small deadzone.

Since stacked bilayers are notoriously hard to label correctly in other methods, we designed a test system containing two closely stacked and intercalating bilayers of POPC (on average two hydration layers of CG water beads, Fig. 3 E). The close proximity and or intercalation of the bilayers did not hinder correct leaflet segmentation at all over a trajectory spanning 10 µs (Fig. 3 F).

The final pure bilayer system tested is a bilayer in which we artificially opened and closed a pore using a biasing potential. We used the leaflet segmentation to detect the instances of pore opening and closing (Fig. 3 G and H, SI video 1).

Consistent leaflet segmentation was successfully obtained using the default settings for all pure bilayer systems even in the presence of cholesterol or intercalating headgroups. The leaflet segmentation can also be used for flip-flop analysis and the detection of a toroidal pore with high accuracy.

Complex bilayers with proteins

After confirming that our leaflet segmentation performs well on simple lipid-only bilayer systems, we went one step further and considered systems of higher compositional complexity.

As for the lipid only segmentation, segmentation of a simple bilayer containing POPC lipids and a single protein is perfect using the default settings (Fig. 4 A, SI video 2). Next, we tested a crowded plasma membrane model. This model contains a high variety of lipids (∼ 60) and a wide range of proteins in a planar membrane. [START_REF] Ingolfsson | Lipid Organization of the Plasma Membrane[END_REF] Segmentation resulted in a nearly perfect assignment (Fig. 4 B). However there are 12 lipids which remain unlabeled (Fig. 4 B zoom box), caused by an exclusion region. Setting the force-segmentation cutoff radius slightly higher resolved the issue, but the current result is shown to demonstrate the quality of the default settings. Missing 12 lipids out of tens of thousands of lipids would probably not affect the analysis to any serious degree. However, if perfection is required some tweaking might be needed. For example the force-segmentation range or the binning distance could be altered as well as the minimum cluster size.

Based on these results we conclude that MDVoxelSegmentation works well on systems with a complex composition even in the presence of protein(s).

Curved and tight geometry

Up to now all our test systems were relatively flat, therefore as our final test set, we introduce heavy curvature as well as closely packed segmentation. In our previous work we showed the phase transition of a stacked bilayer system to an inverted hexagonal phase (H II ), [START_REF] Bruininks | A molecular view on the escape of lipoplexed DNA from the endosome[END_REF][START_REF] Bruininks | Methods in Molecular Biology[END_REF] which was used as a system which should push the time consistency part of our algorithm, as well as the agonistic behaviour with respect to the amount of segments. A snapshot of the first intermediate and last frame of the trajectory is shown in Fig. 5 A. The corresponding graph

shows that the segmentation is extremely stable over time. A link to the video of the full phase transition is given in the SI as video 3. During transition the following happens: First the leaflets which share a common water reservoir (Fig. 5 A I) fuse together by multiple stalk formation (Fig. 5 A II). The leaflets are then split into smaller segments, more closely resembling the classic view of the H II phase. However, due to the high concentration of water in this system not all aquatic channels are singular and some remain connected through water pores perpendicular to the view normal (Fig. 5 A III). The final result is a stable H II phase of channels with a varying degree of interconnection. Note that the visualization of the segments (which is part of the output, as can be seen in the Method section) immediately makes it clear which channels are connected and which not.

In our previous work we showed how transfection with lipoplexes can be simulated using the Martini force field. [START_REF] Bruininks | A molecular view on the escape of lipoplexed DNA from the endosome[END_REF] Such lipoplexes are complexes of lipids and dsDNA. The dsDNA resides inside the inverted hexagonal phase of the lipoplex lipids. Upon fusion with a membrane the dsDNA is released into the cytosol. The tight packing of the lipids with the dsDNA in combination with the change in topology during transfection, makes this another challenging system to analyse. A render of the leaflet segmentation is shown in Fig. 5 B. A link to the video of the complete transfection process is available in the SI video 4. Although the segmentation is pretty good, it does contain some noise, illustrated by minor flickering in the video. This noise could be removed with a noise filter, but we left it in for fair comparison and to create a realistic reference for the quality of the output under comparable conditions. As for pores in the bilayer, we can use the amount of leaflet segments identified in the simulation as a reliable metric for fusion events (Fig. 5 B, graph). We did find that decreasing the resolution of the leaflet segmentation for complex geometry resulted in a steep decrease in quality of segmentation. For more simple geometry 1 nm voxels without hyper-resolution resulted in acceptable segmentation. However, for these systems 0.5 nm voxels with hyper-resolution is required.

The presence of (high) curvature does not hinder correct leaflet segmentation, even if the lipid arrangement does not follow a strict bilayer definition or if the amount of segments is highly dynamic.

Further tests with complex and real world systems

The leaflet segmentation was tested on three additional complex systems to obtain familiarity with the quality. The first system contains a small fatty acid vesicle (Fig. 6, SI video 5).

Around the vesicle fatty acids are suspended in a water buffer. The goal was to measure the imbalance of molecules in the inner and outer leaflet of the vesicle as it absorbs the surrounding fatty acids in solution. As the outer leaflet of the vesicle merges with more and more small micelles in solution the imbalance increases, potentially resulting in division.

Some flip-flopping occurs to release the tension due to the imbalance between the inner and outer leaflet. However, there are two moments at which the inner and outer leaflet of the vesicle temporarily fuse by means of a toroidal pore. These pores quickly cause a significant amount of acyl chains to migrate from the outer to the inner leaflet, indicated by a jump in the inner leaflet's size (Fig. 6 A). Passive flip-flop occurs all the time and is indicated by the more continuous increase in size of the inner leaflet. This example furthermore illustrates the use of our segmentation tool to track flip-flop and pore formation events with high resolution.

The second test was to assign the leaflets in the showcase mitochondrion model presented in the work by Pezeshkian et al. [START_REF] Pezeshkian | Backmapping triangulated surfaces to coarse-grained membrane models[END_REF] For the mitochondrion, we expected to observe four separated leaflets, namely the separate leaflets of both the inner and outer membrane.

However, after leaflet segmentation of the complete mitochondrion we only observed three segments. To further investigate this issue we cut the mitochondria in twenty slices. Each piece was segmented individually (Fig. 6 B). This revealed that there are issues in two of the twenty slices. When these two problematic slices are excluded the expected four leaflets were found by segmentation analysis. This indicates unexpected overlap is occurring in those two slices, likely as a result of stalk formation between closely apposed membrane segments.

Further subdivision could be utilized to pinpoint the exact location, although doing so by hand becomes rather cumbersome. Nevertheless this example shows that large structures can be successfully handled by our algorithm and the results can be used for quality control of systems too large for visual inspection by traditional means.

Finally, in our previous publication on lipoplexes we described the existence of connective channels in an overly hydrated inverted hexagonal lipid phase, shown in the example above on formation of an H II phase(Fig. 5 A). [START_REF] Bruininks | A molecular view on the escape of lipoplexed DNA from the endosome[END_REF] These channels probably have huge effects on the fusion kinetics of such phases. Related work was performed at atomistic scale by Ramezanpour et al., in which they did not describe any of such channels. [START_REF] Ramezanpour | Structural Properties of Inverted Hexagonal Phase: A Hybrid Computational and Experimental Approach[END_REF] After contacting the group and presenting this paradox to them, we were kindly allowed access to their data. We used MDVoxelSegmentation to find any possible unexpected and undescribed connective channels.

The channels were confirmed to also exist in the atomistic simulations, mostly forming at the mid to end of the trajectory of their lipid system with the highest hydration (Fig. 6 C). The channel was stable, but difficult to spot using the conventional tools available. This due to the fact that it spanned over the PBC, a common problem in periodic systems.

MDVoxelSegmentation was able to classify the channel, as it supports all PBC used for molecular simulations. Checking the marked frames by eye confirmed that the results of MDVoxelSegmentation were correct, the authors agreed after reevaluating their data. The formation of these pores was not observed at lower temperatures, suggesting that there is a kinetic factor and that their simulations might benefit from extension.

Together, these examples show that MDVoxelSegmentation is a fast and useful tool for spotting complex segmentation properties, otherwise easily missed or extremely cumbersome to analyse.

Performance

To get some insight in the performance of the current MDVoxelSegmentation tool we timed the segmentation for some of the test systems using a varying amount of threads. The amount of threads were varied to test the parallel processing of MDVoxelSegmentation. In all these tests we turned off force-segmentation, for it can have a large impact on performance based on the exact system (especially at higher maximum recursion depth). The results are summarized in Fig. 7. From the single thread values we can deduce that our algorithm Four distinct leaflets can be traced (green, yellow, purple, red) if the intersection zones are not included in segmentation (white). If the white zones are included in the segmentation, the pink and red leaflets are fused together. (C) Snapshots of the atomistic POPC system in the H II phase as previously published. [START_REF] Ramezanpour | Structural Properties of Inverted Hexagonal Phase: A Hybrid Computational and Experimental Approach[END_REF] Segmentation analysis reveals previously unnoticed connections between the initial four channels at the end of the simulation. The width of a segment in the graph indicates its relative size.

indeed scales pretty much linearly with some variance due to packing and complexity (table 1). Systems containing thousands to millions of particles of interest will roughly take between less than seconds to minutes per frame respectively. The biggest system we timed was the 83 million particle model of a mitochondrion 10, which took over 4 hours for a single frame.

The system of the mitochondrion was the only system which could not be handled with 16GB of RAM and needed a specialized node to run. Its memory consumption peaked at 150GB.

Figure 7: Single frame scalability. The time it takes to segment a single frame with increasing particle count. Exact data is represented in table 1. Some variation from linearity is expected due to the fact that the specific geometry of a system can have an impact on the complexity.

Discussion

Future development of MDVoxelSegmentation

In this manuscript we show the use of voxel based approaches for analysis of lipid and surfactant based systems. However, our implementation is not optimal yet. What we hope to [START_REF] Humphrey | Visual molecular dynamics[END_REF] Pymol, [START_REF] Llc | The PyMOL Molecular Graphics System[END_REF] UnityMol 24 ) to embrace voxels and add user intuitive voxel APIs with common components analysis to their packages. VMD already has a well optimized grid density function, which could be made more accessible for users, or even contain some GUI. The VMD code could probably perform all operations in this manuscript at runtime for systems of considerable size due to their efficient usage of GPU(s). How the segmentation API could best be merged into the current selection syntax is not trivial, and will require careful consideration.

Potential extra features

As shown here, leaflet segmentation can be solved with high fidelity, using only local checks, and chronological segmentation. We showed accurate segmentation of a wide range of lipid systems with or without the presence of cholesterol, protein and or (extreme) curvature.

Although the focus of this manuscript was lipid leaflet detection, we would like to point at some possible additional features that further broaden the scope of the applications. One of those features is the tracking of membrane pores and stalks. This should use the same basic voxel operations such as erosion, growth and common neighbor segmentation, but results in tracking abstract segments such as pores. This would allow us to treat pores as segments, propagating them through the trajectory using the same selection syntax and data structures.

Current attempts show to be promising, although this pore segmentation appears to rely on a good understanding of subdividing the segmentation results. Some preliminary attempts of the porefinding can be found in the MDVoxelSegmentation github pore branches.

Another possible application is in the analysis of interfaces and the manner in which they percolate between set targets. In this case, the interphase itself might be the entity tracked over time using a distance query on multiple segments.

Detecting lipid lateral phase separation is another problem which we think could be potentially tackled with our tool. In this case we should not map lipid densities to voxels, but lipid compositions. In the voxeled compositions a 3D edge detection algorithm could be used to find transitions in local composition. Since we know that a phase its composition should be constant in the bulk region, we can start with finding volumes which are constant in composition by taking the derivative. This would work similar to the work of Sodt et al.

where they demonstrate such approach on a 2D grid. [START_REF] Sodt | The Molecular Structure of the Liquid-Ordered Phase of Lipid Bilayers[END_REF] Probably some form of compositional smoothing or dimensionality reduction would be required to reduce noise in the spatial composition of complex systems.

Finally, the analysis is generic enough to be used practically without adaptations on lipid related systems such as amphipathic proteins or fatty acids, as seen in the SI videos 5 and 6, thus showing promise to be used as a more general procedure for amphipathic systems.

To conclude

We showed that consistent and accurate leaflet segmentation can be obtained with our sequential segmentation algorithm. The algorithm proved to be successful on a wide range of lipid and lipid related systems. These systems contained both high compositional and topological complexity which were not segmentable with any of the currently available methods known to us. To achieve this we used a voxel data structure, which allows for fast local queries, segmenting each frame individually. To achieve consistent segmentation over multiple frames, this was followed by set theory, making use of the Jaccard index. In our tests we

show many examples of CG Martini lipid systems, mostly performed with the same parameter settings. These settings show very few artefacts indicating that in many cases no time needs to be spent on parameter tweaking. Besides CG systems we also looked at an atomistic system for which segmentation showed to be of high quality with only minor changes needed to the default settings.

Our algorithm practically scales linearly with the increase in particles and is able to handle millions of particles in minutes on a desktop machine. However, performance gains of orders of magnitude could probably be achieved by proper SIMD/GPU implementation [START_REF] Cabaret | A review of world's fastest connected component labeling algorithms: Speed and energy estimation[END_REF] and careful memory management using more advanced voxel data structures16-18. The resulting segmentation is easy to work with in python or VMD and can directly be used for visualization due to its supported segmentation data format.

The MDVoxelSegmentation package will be expanded by our own attempts to tackle problems within this framework. We do invite interested researchers to contribute to the further development of the tool. To add your own segmentation routines to the package, simply open a ticket at our github and indicate your interest.

(https://github.com/marrink-lab/MDVoxelSegmentation)

Methods

Molecular dynamics simulations

To simulate the lipid only systems we made use of GROMACS 2019.1. [START_REF] Abraham | GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers[END_REF] Martini 2 was used as the force field. [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF] All lipid parameters were taken from cgmartini.nl and initial configurations were generated using insane. [START_REF] Wassenaar | Computational Lipidomics with insane: A Versatile Tool for Generating Custom Membranes for Molecular Simulations[END_REF] Each simulation contained 0.15 M NaCl and 10% antifreeze (WF) to prevent unnatural freezing of the water as described in the original Martini paper. [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF] Rectangular periodic boundary conditions were employed to prevent artificial wall effects (xy 10 nm). The steepest descent algorithm was used for energy minimization (1,000 steps) and equilibration was performed using the default Martini settings with a time step of 2 fs for 25,000 steps. [START_REF] Marrink | The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations[END_REF][START_REF] De Jong | Martini straight: Boosting performance using a shorter cutoff and GPUs[END_REF] The Verlet cut-off scheme was used with a 1.1 cut-off for both the coulombic (reaction-field) and van der Waals interactions. The v-rescale scheme (tau-p=1 ps) was used for the thermostat at 310 K, coupling the lipids and solvent/ions in two separate groups. Pressure coupling during equilibration was performed using the Berendsen barostat for semi-anisotropic systems (tau-p=3 ps, compressibility=3e-4 bar-1, ref-p=1.0 bar in each dimension). [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] The production run made use of a 20 fs time step and the pressure coupling was switched to Parrinello-Rahman (tau-p=12). [START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] To simulate pores in the DLPC bilayer, positional restraints were used in the form of a cylindrical biasing potential acting on the C1A and C1B tail beads (GROMACS potential function 2, shape 8, hole radius -1.2 nm, force constant 1,000 kJ mol-1). The negative sign in front of the hole radius indicates an inversion of the biasing potential in GROMACS.

A reference file was created with all beads at the same position, acting as the pore center reference position (-r flag in GROMACS). The pore formation was simulated by running a production run for 10 ns without the biasing potential, then the potential was turned on for another 10 ns and turned off again for the final 10 ns.

Leaflet segmentation

We used the same default settings for all the presented segmentation unless specified differently. These settings used a binning resolution of 5 Å and made use of hyper-resolution of 0.5 (adding points half the resolution away from the source point). The headgroups, linkers, tails and exclusions were as specified in the default selection input. The minimum size of a segment was 50 particles and iterative force-segmentation was turned on with a maximum distance of 20 Å. An identity threshold of 0.618 was used for the Jaccard threshold. All keyword arguments represent default values which do not need to be specified except for the reference frame and trajectory file (any MDAnalysis supported MD format). For the segmentation of the atomistic H II system hyper-resolution was turned off (force segmentation 0). In general hyper-resolution is not required for atomistic systems.

The respective headgroups, tails and linkers selections were specified in the selections input file:

[charmm_heads]

"name N P C12 C11 O11 O12 O13 O14" Flip-flop comparison with smaller deadzone The role of the Jaccard threshold

For temporal segmentation a Jaccard threshold had to be set. This threshold defines when two sets are similar enough to be a candidate for shared identity. Picking a very low threshold results in an overzealous attempt to never generate a new set, where a very high threshold results in a flood of new identities -a single lipid flip-flopping can cause a loss of identity.

We found that values around the golden ratio correspond well with our human intuition.

Figure 10: Comparing Jaccard thresholds. Three values for the Jaccard threshold (0.000, 0.618, 0.999) were tested on two bilayer systems. A bilayer with a toroidal pore (A). A bilayer with cholesterol (B). The Jaccard threshold of 0.618 gives the desired results in both systems.

Setting the threshold to 0 results in incorrect temporal segmentation of system A, although system B is fine. Setting the threshold to 0.999 results in bad temporal segmentation in both systems. Segment IDs beyond 50 were truncated.

Supporting Videos

Supporting Video 1: Bilayer pore bilayer pore The opening and closing of an artificial pore in a DLPC bilayer. Right before the stable toroidal pore is formed, a short lived water pore spans the bilayer. This pore was visualized using the following selection syntax in VMD: '(name W) and (within 20 of name PO4 and user 1) and (within 20 name PO4 and user 2)'. User 1 and 2 are the respective leaflets. We can visualize water pores connecting two leaflets and segment them. However, can happen to restore the balance. The acyl chains flip-flop from the outside to the inside.

The leaflets temporarily fuse to balance the amount of acyl chains in the inner and outer leaflet, or the disbalance causes the vesicle to divide. MDVoxelSegmentation was used to analyse the imbalance as well as to monitor the feeding rate. Both flip-flopping and toroidal pore formation takes place. No divisions occur.

Supporting Video 6: Segmentation of other amphipathic systems amphipathic peptides With trivial adaptations of the input GRO file we used our method to perform segmentation of amphipathic peptides. Due to their amphipathic nature, micelles from such peptides can touch but not fuse (like lipid leaflets). MDVoxelSegmentation shows to be suited to be used for such systems.

Figure 1 :

 1 Figure 1: Spatial leaflet segmentation. A 2D schematic representation of the leaflet segmentation of two stacked bilayers in a periodic box. The upper bilayer contains a protein channel and the lower bilayer contains a diving lipid. The protein is used as exclusion selection.

Figure 2 :

 2 Figure 2: Finding the identity of the segments in the new frame. (A) The ordered 1 segment has a best match with unordered 2. (B) The ordered 1 and 2 are merged into one segment, with ordered 2 providing the identity for unordered 1. (C) The ordered 2 is split into two new labels, unordered 1 and 2. Ordered 0 is always matched as the new segment zero, as this is the segment which contains all non-segmented particles. Capital N indicates the total number of segments in each frame.

Figure 3 :

 3 Figure 3: Leaflet segmentation and analysis for pure bilayer systems. (A, B) Snapshots of a DOPC/CHOL(4:1) bilayer with and without force-segmentation. Some lipids were unassigned without force-segmentation due to diving/occlusion (A, yellow). Using forcesegmentation non-segmented lipids are assigned to a segment as a post-processing step (B). (C) The flip-flop traces of a single cholesterol in the DOPC/CHOL bilayer, comparing the z-height and segmentation flip-flop approaches. Whenever the orange line is not visible it lies perfectly behind the green line. (D) A comparison of the up-and downward cholesterol flipflop events between the z-height and segmentation based flip-flop analyses. A more intense single dot represents more than one cholesterol flip-flopping at that given instant. (E) A snapshot of closely stacked and intercalating bilayers of DLPC. (F) Leaflet segmentation of the intercalating bilayers over 10 µs. The width of a segment represents its relative size (G) A snapshot of a DLPC bilayer with a toroidal pore (frame 750 of H). (H) Leaflet segmentation of pore formation. The pore forming biasing potential was turned on between frame 500 and 1000. The width of a segment represents its relative size.

Figure 4 :

 4 Figure 4: Leaflet segmentation in the presence of protein. (A) Adding protein (lime) to a bilayer does not interfere with the leaflet assignment of a CG POPC bilayer (red/blue).A video of the bilayer protein system can be found in the SI. (B) The bottom and top leaflet of the plasma membrane. Even extremely complex bilayers such as the plasma membrane including many different proteins (transparent/lime) and lipids (red/blue) do not hinder leaflet assignment. However, there is a small cluster of non-assigned lipids (yellow) around one of the proteins (zoom box).

Figure 5 :

 5 Figure 5: Leaflet segmentation and analysis for both tight and curved complex geometry. (A) The phase transition from a multilamellar bilayer system to an H II phase. Snapshots of the first, intermediate and last segmented frame of a 1 µs trajectory. The accessory video is available in the SI at supporting video 3. The corresponding graph shows all segments formed in chronological order with the snapshots indicated with the dotted lines. The width of a segment indicates its relative size. (B) Snapshots of a lipoplex transfecting dsDNA over a bilayer (first, intermediate and last segmented frame). The inner core (orange, yellow, green, gray) is an intact H II phase containing the dsDNA (not visible), whereas the upper bilayer leaflet (red) has already fused to the outer leaflet of the lipoplex by means of a fusion stalk (I). During transfection the inner core segments fuse with the bottom leaflet of the bilayer (II), resulting in a single flat bilayer with all the dsDNA on the other side (III). The transfection can be tracked in detail using the segmentation graph. The width of a segment indicates its relative size. 16

Figure 6 :

 6 Figure 6: Additional systems. (A) A snapshot of a fatty acid vesicle (red/grey), surrounded by individual fatty acids in the medium (green). The amount of fatty acids in the inner and outer leaflet of the vesicle over time (graph). A video of the fatty acid dinner is available in SI supporting video 5. (B) A snapshot of the segmented mitochondrion (˜83 million beads).Four distinct leaflets can be traced (green, yellow, purple, red) if the intersection zones are not included in segmentation (white). If the white zones are included in the segmentation, the pink and red leaflets are fused together. (C) Snapshots of the atomistic POPC system in the H II phase as previously published.[START_REF] Ramezanpour | Structural Properties of Inverted Hexagonal Phase: A Hybrid Computational and Experimental Approach[END_REF] Segmentation analysis reveals previously unnoticed connections between the initial four channels at the end of the simulation. The width of a segment in the graph indicates its relative size.

5 --recursion_depth 10 --force_segmentation 20 -

 51020 mdvseg -f start.gro -x trajectory.xtc --headgroups martini_heads --linkergroup martini_linkers --tailgroups martini_tails --exclusiongroup martini_proteins --resolution 5 --hyper_resolution 0.-force_information False --minimum_size 50 --begin 0 --end None --stride 1 --threads 12 --bit_size 32 --output clusters --verbose False

  C23 C24 C25 C26 C27 C28 C29 C210 C211 C212 C213 C214 C215 C216 C217 C218 C32 C33 C34 C35 C36 C37 C38 C39 C310 C311 C312 C313 C314 C315 C316 C317 C318"The AA terminal input was: mdvseg -f start.gro -x trajectory.xtc -hg charmm_heads -lg charmm_linkers -tg charmm_tails -hres 0

Figure 8 :

 8 Figure 8: Qualitative comparison of leaflet segmentation algorithms. A variety of lipid systems was segmented using the default options of MDAleaflets, FATSLiM and MDVox-elSegmentation (MDVseg). Solid images indicate reliable output. Faded images represent the best case output for unreliable segmentation. The headgroups were selected using the default MDVoxelSegmentation input selection. The large stop sign indicates no output was generated at all during segmentation. Unreliable segmentation is indicated with the fingers crossed symbol. *The speed is indicated relative to the MDVoxelSegmentation performance. **Segmentation can be extremely sensitive to diving lipids. ***The option for sparse matrices does exist, however, it is not the default behavior and the calculation speed suffered severely. 28

Figure 9 :

 9 Figure9: Comparing flip-flop detection with a small deadzone. Reducing the deadzone for the z-height based flip-flop registration to 0.3 showed that the flip-flops detected by the MDVoxelSegmentation approach is a near perfect subset of the z-height approach with a small deadzone. Although the exact moment of the flip-flop is sometimes shifted by a few frames between the two approaches.

Table 1 :

 1 Performance of MDVoxelSegmentation. All systems were run on the same Ubuntu desktop except for the mitochondrion. The test machine contained a Ryzen 1600 (6 cores, hyper threading 3.8 GHz) CPU and 16GB of DDR4 memory (2x 8GB 2.133 MHz). No GPU acceleration is currently implemented. The MDVoxelSegmentation parameters were the same as stated in the examples above except for force-segmentation, which was turned off (-fs 0). Timing was performed using the default Ubuntu time functionality. The CPU for the mitochondrion system ran at 2.4 GHz for a single thread.

	System	Particles Frames Threads Real time User time Sec/Frame
	Bilayer pore	10,436	1,500	12	2m25s	21m41s	0.096
	Fatty acid growth	66,278	1,000	12	3m35s	27m37s	0.215
	Bilayer pore	10,436	1,500	1	13m13s	15m23s	0.529
	Fatty acid growth	66,278	1,000	1	16m20s	18m30s	0.980
	Plasma membrane 1,343,450	1	1	4m35s	4m36s	275
	Mitochondrion	83,288,300	1	1	260m18s	-	15,618
	achieve is to stimulate developers of MD software which have a large user base (MDAnal-
	ysis, 12 VMD,

Acknowledgement

Bart M. H. Bruininks thanks Manuel N. Melo, Jonathan Barnoud, Weria Pezeshkian, John Stone and the rest of the VMD developers, as well as Elmar Eisemann for the many fruitful discussions. We also thank Melanie König (fatty acids), Haleh Abdizadeh (plasma membrane), Josef Melcr (amphipathic peptides), Weria Pezeshkian (mitochondrion) and Peter

Lipid flip-flop

The average z-height of the cholesterol headgroups was subtracted from the headgroup zposition. This resulted in headgroups lying either above or below zero. Headgroups with a normalized z-position less than 0.6 nm away from 0 were labeled either up or down (discretization of the z-height). Differentiating the individual lipid traces results in peaks at flip-flop events. These peaks were counted, registering if the flip occurred up-or downwards. This calculation was repeated for all cholesterols in all frames and a graph was made using python3 matplotlib. [START_REF] Hunter | Matplotlib: A 2D Graphics Environment[END_REF] The z-height and flip-flop analysis is available at the MDVoxelSegmentation github.

To show the flip-flopping of cholesterol using our leaflet segmentation we used the default settings except for the force-segmentation which was turned off. Lipids which were not assigned to segment 0 were assigned to their previous segment. A graph was made using matplotlib and the plotting script included with mdvseg. The leaflet flip-flop analysis is available at the MDVoxelSegmentation github.

Visualization

All figures and videos were rendered using VMD and a little TCL script of MDVoxelSegmentation to load the segmentation array. Examples of TCL scripts are provided together with the code at: https://github.com/marrink-lab/MDVoxelSegmentation Supporting Information Available Supporting Figures

Comparing leaflet detection algorithms

A qualitative comparison of segmentation quality between MDAleaflets, 12 FATSLiM 14 and MDVoxelSegmentation was performed to explore the strengths and weaknesses of each approach.

we lose all information of the pore as soon as the leaflets become continuous (as illustrated in the video).

Supporting Video 2: Protein bilayer protein bilayer Stable segmentation of a bilayer in the presence of a protein. Close to the protein some lipids can be observed which are diving relatively deep under the surface, however, this diving is handled without missegmentation.

Supporting Video 3: Lipid phase transition phase transition Lipid segmentation of a phase transition from dehydrated stacked bilayers to inverted hexagonal (left). Water densities (right).

Supporting Video 4: Lipoplex transfection

transfection A lipoplex was placed on top of a model endosomal membrane. The simulation was started in the presence of a fusion stalk, connecting the outer leaflet of the lipoplex coat with the upper leaflet of the endosomal membrane (red, transparent). Over time the internal compartments of the lipoplex (orange, yellow, pink, green) fuse with the bottom endosomal leaflet (red, transparent). At the end of the simulation all dsDNA from inside the lipoplex was transfected over the endosomal membrane. Some flickering is observed, but overall the segmentation is very stable and flickering could easily be removed with a noise filter.

Supporting Video 5: Acyl chain dinner acyl chains, acyl chains dof In an attempt to understand early life 'cell' division we placed a small vesicle (acyl chains) in the presence of food (also acyl chains). The vesicle adds the food from its surroundings to its own composition, but without internal reorganization the food is only added to the outer leaflet. With a large enough leaflet disbalance three things

Graphical TOC Entry

Some journals require a graphical entry for the Table of Contents. This should be laid out "print ready" so that the sizing of the text is correct. Inside the tocentry environment, the font used is Helvetica 8 pt, as required by Journal of the American Chemical Society. The surrounding frame is 9 cm by 3.5 cm, which is the maximum permitted for Journal of the American Chemical Society graphical table of content entries. The box will not resize if the content is too big: instead it will overflow the edge of the box. This box and the associated title will always be printed on a separate page at the end of the document.