{'fraction_non_alphanumeric': 0.025533163764429663, 'fraction_numerical': 0.009261070892845497, 'mean_word_length': 4.769896519285043, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Efficient Algorithm for Topological Characterisation of Worm-Like and Branched Micelle Structures from Simulations

Many surfactant-based formulations are utilised in industry as they produce desirable visco-elastic properties at low-concentrations. These properties are due to the presence of worm-like micelles (WLM) and, as a result, understanding the processes that lead to WLM formation is of significant interest. Various experimental techniques have been applied with some success to this problem but can encounter issues probing key microscopic characteristics or the specific regimes of interest. The complementary use of computer simulations could provide an alternate route to accessing their structural and dynamic behaviour. However, few computational methods exist for measuring key characteristics of WLMs formed in particle simulations. Further, their mathematical formulation are challenged by WLMs with sharp curvature profiles or density fluctuations along the backbone. Here we present a new topological algorithm for identifying and characterising WLMs micelles in particle simulations which has desirable mathematical properties that address short-comings in previous techniques. We apply the algorithm to the case of Sodium dodecyl sulfate (SDS) micelles to demonstrate how it can be used to construct a comprehensive topological characterisation of the observed structures.

lead to WLM formation is of significant interest. Various experimental techniques have been applied with some success to this problem but can encounter issues probing key microscopic characteristics or the specific regimes of interest. The complementary use of computer simulations could provide an alternate route to accessing their structural and dynamic behaviour. However, few computational methods exist for measuring key characteristics of WLMs formed in particle simulations. Further, their mathematical formulation are challenged by WLMs with sharp curvature profiles or density fluctuations along the backbone. Here we present a new topological algorithm for identifying and characterising WLMs micelles in particle simulations which has desirable mathematical properties that address short-comings in previous techniques. We apply the algorithm to the case of Sodium dodecyl sulfate (SDS) micelles to demonstrate how it can be used to construct a comprehensive topological characterisation of the observed structures.

Introduction
Micelles play an important role in many processes and products. For example, micelles can be used to control drug release, [1][2][3][4][5][6][7] are key in household cleaning products, 8,9 and can act as friction modifiers in vehicle engines. 10,11 They display a wide variety of fluctuating behaviours and morphologies depending upon the surfactant from which they are assembled and the environment in which they reside. The different structures formed by micelles give rise to multiple possible behaviours including interesting optical and rheological properties.
For example, rod-and worm-like micelles (WLMs) are associated with a significant increase in overall viscosity of a liquid when compared to spherical micelles. WLMs are long flexible cylinders with diameters (for typical surfactants molecules) of approx 3-5 nm and persistence lengths of approx. 20 nm below which they should be considered as rigid rods. These micelles are able to grow to a length of several micrometers at relatively low surfactant concentrations resulting in the occurrence of entanglements and ultimately enhanced visco-elasticity. [12][13][14][15] Under certain conditions WLMs may branch between each other to form interconnected networks. For example, electrostatic screening by counter-ions is known to promote WLM branch formation. [16][17][18] Whilst long entangled micelles lead to increased zero-shear viscosity, branching and network formation is generally believed to lower it. 19,20 Thermal fluctuations can cause micelle branches to slide along the main micelle backbone chain thereby providing an accessible mode of stress relaxation. 19,21 It is believed this in turn leads to a reduction in viscosity when compared to an entangled matrix of WLMs. 17,[21][22][23][24] The formation of WLM structures is dependent upon the underlying structural characteristics of the surfactants molecules from which they are constructed. Surfactants that form WLMs preferentially pack into cylindrical shapes. Additives can significantly alter this preference through, for example, hydrophobic binding and electrostatic screening, leading to changes in WLM structure and rheological properties.
Exploring the behaviour of WLMs can be challenging using experimental methods. Cryo-TEM gives the researcher a window into the micelle structure but this approach can be difficult because of the high viscosities of concentrated surfactant solutions. Pulsed gradient spinecho nuclear magnetic resonance (PGSE-NMR) 25 and multiple scattering approaches 26,27 can supply information on the branching behaviour of micelles. However, branching is most often inferred from rheological observations (i.e., a decrease in the zero shear viscosity).
Computer simulations offer a complementary approach to existing experimental techniques for studying WLMs and are beginning to shine some light onto the structural and dynamic behaviour of WLMs. Self-assembly of surfactant molecules into micelles is difficult to simulate with all-atom methods, such as Molecular Dynamics (MD), due to the high computational cost associated with obtaining appropriate length and time-scales. Consequently only a few micelles can currently be studied with these techniques. As such, coarsegrained (CG) approaches are often applied in the study of micelles, such as the MARTINI force-field. 9,16,18,[28][29][30][31][32][33][34][35][36][37][38][39][40] Dissipative Particle Dynamics (DPD) is another CG approach which is becoming increasingly popular for the simulation of micelle behaviour as it allows signif-icantly longer length and time-scales to be accessed at reduced computational cost (albeit with a potential loss in accuracy). 41-43 DPD employs very soft conservative potentials which facilitates the use of large time-steps. Further details of the DPD method can be found in the literature. [41][42][43][44] A recent review on simulating surfactants and micelle systems (including WLMs) using a variety of particle-based simulation techniques has been presented by Taddese et al. 45 In addition, a review of the application of computer simulation to the study of WLMs can be found in reference. 46 A number of rudimentary analytical approaches have been used to characterise micelles including aggregation numbers and size distributions, radius of gyration, asphericity 47 which are often adequate for small micelles with simple topological structures but less appropriate for larger and more topologically complex structures such as WLMs. Here quantification of branch points, micelle backbone length and end-cap formation are more important but far less methods exist in the literature for capturing these properties. A notable exception is the work of Dhakal and Sureshkumar who have employed a contour mesh measurement approach to calculate a micelle backbone length scale and identify micelle branch points and end-caps. 16 This model was later successfully employed to examine the stress relaxation modes of branched micelles. 21 Nevertheless, their contour length measurement algorithm is mathematically limited to piecewise linear backbone length metrics where each constituent line is of equal length. A more comprehensive topological framework is required to accurately represent the geometric complexity of micelles with sharp curvature profiles or density fluctuations along the backbone. In particular, these latter considerations cannot be overlooked when one wishes to simulate multi-component industrial formulations where surfactant polydispersity, and the presence of co-surfactants and modifiers potentially require additional scrutiny.
The remainder of this article is arranged as follows. First, we give an overview of the topological algorithm applied to characterise the salient geometric length scales of WLMs and branches. Next, we detail the DPD simulation protocol and SDS micelle system test case considered in this study. We then present our observations for the SDS test case. The paper concludes with a discussion on the developed topological algorithm including how some of the computed metrics can be utilised as input parameters to macroscopic micelle viscoelasticity models.

Algorithmic Design
The framework behind the topological characterisation algorithm introduced in this article is outlined in Fig. 1. First, individual micelles are identified using trajectory data from DPD simulations (See Section 3.1). Each micelle can be represented in the form of a point cloud (and neighbourhood graph) which is the set of 3D spatial coordinates of each molecular bead in the cluster (Fig. 1a). Next, the point cloud structure undergoes an operation known as Fiedler vector segmentation which partitions the structure into a number of distinct regions with corresponding centre points (Fig. 1b). This exercise is repeated until an optimal spectral clustering is achieved which captures important features such as branch points and loops without overfitting (Fig. 2). Finally, these points are joined together to produce a piecewise linear (PL) skeleton of the micelle (Fig. 1c) composed of a set of consecutive linear segments of varying individual lengths. This key step allows us to calculate a series of structural micelle characteristics such as the micelle backbone length, the number of branch points and the Euclidean distance between end-caps.
A description of the theory underpinning the mathematical framework is presented in the Supporting Information.

Spectral Splits
The spectral gap, denoted here by λ, is equal to the first eigenvalue of the discrete Laplacian. 48 The general protocol is to recursively split sets at their spectral gaps (with further information outlined in the Supporting Information). This process works well for large WLMs; however, overfitting from unnecessary cuts leads to excess length and inaccurate branch points in smaller micelles without clear WLM structure. To avoid this, we introduce a bound λ max on the spectral gap λ of each potential split which ensures a sufficiently sparse cut. We note that recursion introduces new eigenvalues for each segment generated and new spectral gaps. On these data sets, we have observed that the spectral gaps increase during recursion.
There are two fundamental properties of a good split that can be revealed by the spectral    Empirically, we have discovered that terminating the recursion when any newly computed spectral gap value exceeds λ max = 0.015 has been a good stopping criterion over the data tested. Since this was determined empirically, there is no claim that the value of 0.015 is optimal.

Segmentation
The recursive procedure Segment is detailed in Algorithm 1. For each recursive step the spectral gap λ is used as a stopping condition so that recursion can only continue when λ > λ max . To avoid overfitting we impose a soft maximum recursion depth k max . To avoid underfitting, we set a lower bound λ min on the spectral gap so that k max only applies when λ min ≤ λ ≤ λ max . It can be shown that the gap of each of the resulting subgraphs must be greater than λ and continued partitioning will converge to λ min .
Input graph G and depth counter k. L ← NormalizedLaplacian(G) Compute the normalized Laplacian of G.
Eigenvalues λ i and vectors w i sorted by The Split procedure, detailed in Algorithm 2, takes a graph G and partitions its vertices by the vector f . For each split the resulting graph partition has vertices corresponding to subgraphs G 1 and G 2 which are stored in a tree structure T with children( If G is connected, a single split will result in one edge added to the graph partition which corresponds to a link subgraph G l that associates its children. This association captures the vertices in the graph partition which are split by G l and will be denoted by nbr(G l ) so that nbr(G l ) = {G 1 , G 2 } initially. The asymptotic complexity is bounded above, quadratically, 49 by the eigendecomposition of a graph Laplacian, but implementation optimizations greatly improve the amortized running time in practice.
Although each step of the segmentation resembles a graph partition the resulting structure is not, strictly speaking, a graph partition. This is due to the fact that when child subgraphs G 1 and G 2 are split further, nbr(G l ) is updated to include the children of G 1 and for all L j ∈ links(T ), where the intersection is taken as the intersection of vertex sets.
Note that a link node may associate more than two leaf nodes, indicating a branch point in the skeleton. An additional step merges intersecting links under certain conditions. The resulting structure is therefore not a graph partition as we treat the links as vertices of the skeleton which may have degree greater than two.

Algorithm 2
The split procedure. Require: Tree T with subgraphs as nodes.
Stores parent / child relationships and links between leaf nodes.
Merge links G l and L. end if end for return G 1 , G 2 end function When the segmentation procedure terminates the leaf nodes of the resulting tree structure correspond to a collection of subgraphs with vertex sets that partition the input vertices.

Skeletonization
Recall that each vertex v i in a micelle graph G = (V, E) is associated with a point p i ∈ R 3 that gives the location of the corresponding DPD bead in the micelle. Therefore, the vertex set of a subgraph G = (V , E ) of G can be associated with a collection of points in R 3 , Let φ be the centroid function, which maps subsets of R 3 to a single point, by where |P | denotes the cardinality of P .
For a micelle graph G = (V, E) the segmentation procedure yields a tree T such that and edge set The skeleton of a segmentation graph by a function φ is the graph and corresponds to a collection of line segments in R 3 drawn between the centroids of points in the leaves and links of T .

Estimating the Full Length
We first produce an approximating PL curve that connects the centroids. This curve remains within the interior of the point cloud, so its length is an underestimation. We add additional cap nodes to approximate the full length. At each terminal point of the curve, extend the containing segment until it exits the point cloud. The cap node is that exit point. Figure 4 depicts a case in which the cap nodes yield a far more accurate calculation of length.

Skeleton Repair
After initial skeletonization, additional segmentation may be required for more complex micelles, in particular for those with branched structure. Fig. 5a depicts a case in which the skeleton does not accurately capture the micelle's structure, resulting in a missed cap node. This is due to a degenerate case in which an end segment has two connections which are too close (Fig. 5b). The final step repairs any degenerate end segments, ensuring that the resulting skeleton accurately reflects both the length and the branching structure of the micelle (Fig. 5c).   5 shows a representative situation, as the blue node has two neighboring nodes, shown in red. The Euclidean distance between the two red nodes is quite small compared to the distance between them, measured along the PL curve. This intuition is captured as a bound on the ratio of the two distances. In particular, we repair a subset of the PL curve when this ratio is less than 0.2. In Fig. 5c the repair is indicated by a new blue node, as a common point for two new subsets of a new PL curve.

Estimating a Cross-sectional Radius
After the skeleton has been stabilized, an estimate of a cross-sectional radius for each WLM can be computed. The role of this radius is to approximate a tube of fixed radius, around the skeleton, so that all points of the WLM lie within the tube. For each cluster, the minimum Euclidean distance from each cluster point to the local linear segment of the skeleton is computed, producing a local radius for each cluster. The final estimate r cs for the crosssectional radius is the minimum over all these radii as a reliable lower bound.

Simulation Protocol
To explore the capabilities of the topological characterisation algorithm we performed nine

Test Case: SDS Micelle Properties in the Presence of NaCl
In this section we examine the application of the topological algorithm to elucidate the structural evolution of SDS micelles as a function of salt concentration. Fig. 6 depicts the various micelle morphologies extracted from the simulations at different levels of NaCl, where each cluster is coloured as per its shape categorisation derived from the UMMAP package. 47 Our topological algorithm partitions the extended UMMAP shape category (orange micelles in Fig. 6) into WLMs and branched micelles, producing the calculated micelle metrics shown Figure 6: Snapshots of the micelle morphologies of each simulation. The caption of each subplot denotes the salt concentration (w/w %) while the colors denote chemically significant changes in shape. 47 in Fig. 7.
In the absence of salt (Fig. 6a), the electrostatic interactions between the sulphate head groups produces an unscreened effective repulsive barrier. This inhibits cluster aggregation and results in the formation of small spherical clusters where the surface area is minimised in order to reduce repulsive head group electrostatic interactions. 18 As salt is added to the system, it screens these electrostatic interactions leading to micelle elongation (Fig. 6b) and a reduction in their number (Fig. 7e).
As salt concentration is incremented further, the disruption of the electrostatic repulsion facilitates the transition from rod-like to WLM micelles (Fig. 6b & c). Eventually, it is screened to such an extent that branching node points can develop (Fig. 6c & d) as measured in Fig. 7d. That is to say that branch points where three micelle segments intersect become energetically possible. This signifies the emergence of branched micelles and coincides with further reductions in the number of micelles (Fig. 7e) and expanding micelle backbone networks (Fig. 7b). This trend continues until we end up with one extensive branched micelle which is progressively developing more node points (Fig. 7f).

Discussion
Although micelle length is known to be a key metric associated with structural changes in rods and WLMs 60 there has been no obvious way to measure the length of WLMs.
For a cylinder, measuring the axis provides an obvious overall length. Analogously, for WLMs we introduce the PL skeleton as a metric comparable to the axis of the cylinder, while still permitting an easy determination of overall length (by the obvious addition of the lengths of all the sub-segments). This skeleton is similar to the medial axis, 61-66 which plays a prominent role in geometric topology algorithms for many engineering applications.
However, the approach here avoids the topological instability problem where small changes in shape cause significant variations in the computed medial axis. 67 The effectiveness of this length value as a metric is corroborated by the chemical simulations reported in the previous section.
Our procedure contains three parameters controlling spectral segmentation which were introduced to avoid observed problem cases of data fitting and whose values were empirically determined. The lower bound on the spectral gap, λ min avoids under-fitting while the upper bound, λ max , and the maximum recursion limit, k max avoid overfitting. In our simulated experiments we found that a balance of these parameters was essential to maintaining good performance. After tuning these parameters we still encountered cases in which spectral segmentation alone failed. An additional repair and merge step was introduced to handle these, which included branch points with high degree and the presence of cycles in the skeleton graph.
The algorithm was designed from a geometric topology perspective with assumptions made about the shape. We assume that the input graph is connected and represents a solid object, with high connectivity in the interior. This assumption is vital to the effectiveness of spectral segmentation as it allows the boundary and branch points to be clearly identified.
While our empirically determined parameter values work well in practice under these shape assumptions, we hypothesize that these thresholds reflect basic structural characteristics of the micelles and could be derived directly from a more careful chemical analysis.
Our topological algorithm can be successfully applied to convex micelle geometries, WLMs and branched structures, however there are limitations to its range of applicability. In particular, surfactant bilayer and vesicle micelles were not investigated in this study and further work would be required to extend the mathematical framework such that it could correctly characterise such point cloud configurations.
A key motivation for our algorithm is to accurately measure structural characteristics of WLMs that could be used in conjunction with viscosity models to obtain further insight into their rheological properties. For example, based on polymer physics concepts Boek and coworkers developed a framework, MESOWORM, which models WLMs as a continuous string of thin rods, similar to our derivation of a piecewise linear micelle backbone. 68-72 From this they derived an expression for the zero-shear viscosity, η 0 , of unentangled WLMs. 69 The longest relaxation time, τ l , is given as: Here b is the average micelle segment length which can be extracted from our algorithm. ξ is the average friction on a segment and can be estimated as outlined in the literature. 69 The final parameter is the effective relaxation time τ . This is a function of τ l (Eq. 2) and τ b , the characteristic breakage timescale, which equals 1/(k m L). Here L is the average micelle length produced from our algorithm and k m is the micelle breakage rate. Although k m is not measured by our algorithm, the fact it can measure end-caps allows us to detect their formation, and thus provides an avenue for obtaining this value from the same simulations.
We note that the practical use of the above equations will require careful consideration of finite-size effects and related calibration of the simulation box size. These effects directly impact measurements of extrinsic properties like micelle length and can influence intrinsic properties like micelle backbone segments per unit volume. In addition the growth of the micelles may be inhibited by the finite number of surfactants initialised in the simulation box. These issues will be investigated in more depth in future studies of surfactant rheology.
Finally the algorithm's ability to detect micellar phase transitions (see Section 4) could be useful from a theoretical and experimental perspective. Qualitatively, the onset of branching makes additional stress relaxation modes available, 21 lowering viscosity and indicating that the salt concentration corresponding to a viscosity maximum has been passed. Importantly, the onset of branching cannot be detected experimentally and our model could act as a preliminary screening tool for narrowing down the position of the salt curve maximum, and also help to delineate the bounds of applicability of theoretical viscosity models which cannot be employed on branched WLMs.

Conclusion
The work presented here tailors techniques from computational geometric topology to address the problem of measuring the length of worm-like micelles (WLMs) resulting in a novel contribution to their analysis. The developed algorithm takes molecular point clouds as input and produce clusters, as determined from the Fiedler vector of the discrete Laplacian, to form piecewise linear (PL) skeletons, and to assign a comprehensive length to each WLM.
The new length metric allows calculation of micelle-curvature, discriminates worms from rods, and extends directly to branched WLMs. In addition, the algorithm provides information on other interesting characteristics of worms such as the number of end-caps and branch-points, along with the micelle cross-sectional radius. Since a key motivation for characterising WLMs is understanding their rheological behaviour we also discuss the theorized connections between the characteristics accessed by our algorithm and viscosity, and show how the algorithm can provide input parameters for well-established mesoscopic WLM viscosity models.
Importantly, the technique can be used to characterise WLMs produced by any particle

Appendix: Coarse-Grained DPD Model
In this work we adopt the Dissipative Particle Dynamics (DPD) simulation approach in which particles (referred to as DPD beads hereafter) interact via soft repulsions and with local pairwise dissipative and random forces which work together as a thermostat. 42 We do not repeat here details of what is now a standard simulation method. Rather, we point the reader to chapter 17 of the textbook by Frenkel and Smit, 44 and the original DPD literature. 41,42,73 An up-to-date perspective on the DPD methodology has been recently presented by Español and Warren. 43 In the DPD simulations carried out in this article we adopt the sodium dodecyl sulphate (SDS) model developed by by Anderson et al. 74,75 This model has been adopted and extended by a number of groups to explore the behaviour of ionic surfactants. [76][77][78][79] In their approach each DPD bead represents between 1-3 "heavy" atoms (namely, oxygen, Full details of the model and adopted parameters are presented in the SI. In this model a standard set of reduced units is adopted in which the beads have unit mass, the system is governed by temperature k B T = 1 (equivalent to 298 K), and the baseline cut-off distance for the short-range soft pairwise repulsion between solvent beads is set as Supporting Information:

Algorithmic Theory
In the topological characterisation procedure we represent a micelle as a graph S1 associated with a collection of points in R 3 . Formally, a graph is an an ordered pair of sets where V is called the vertex set, with vertices, v i , and E is called the edge set with edges, e i,j . Figure S1 shows a planar graph, with integers indicating vertices and line segments for edges. We consider undirected graphs, so e i,j = e j,i . The degree d i of a vertex v i is the number of edges in E containing v i . The degree matrix of a graph with n vertices is the n × n diagonal matrix D such that The adjacency matrix of G is the n × n symmetric matrix A with entries The graph Laplacian S1-S4 of a graph G is the n × n matrix L = D − A.

S-2
In Figure S1 the vertex labeled with 1 has 2 edges and the degree of vertex 1 is 2, for A standard refinement of the Laplacian defined above which is often used to compare matrices that are of varying size is the normalized symmetric Laplacian, S5 denoted as L sym .
The adjustments are given by entries Throughout, we will denote the normalized laplacian L sym by L.
Let λ 0 ≤ λ 1 ≤ . . . ≤ λ n−1 denote the eigenvalues of L and w i ∈ R n denote the eigenvector corresponding to λ i : that is, Lv i = λ i w i . For a connected graph λ 0 = 0 and λ 1 is the spectral gap which gives a measure of the connectedness of the graph. The eigenvector corresponding to the smallest non-zero eigenvalue, w 1 if the graph is connected, is known as the Fiedler vector. S2,S6,S7 We will assume G is connected and denote f = w 1 and λ = λ 1 to simplify notation.
Vertices will be denoted by v i for i = 1, . . . , n. Entries of f are denoted f i , for i = 1, . . . , n.
The graph G = (V, E) can be split into subgraphs G 1 and G 2 by partitioning V into S-3 and taking G 1 and G 2 as induced subgraphs of G.
Formally, a graph partition is the reduction of a graph G = (V, E) into a collection of smaller graphs with H = (U, L) formed by partitioning its vertices V into mutually exclusive, non-empty, vertex sets U = V i , . . . , V m . That is, each vertex of H is a subgraph H i = (V i , E i ) such that V i ∩ V j = ∅ for all i = j and i=1,...,n V i = V with edges E i consisting of edges in G with both ends in V i . The edges L of the graph partition are given by collections of edges in the initial graph that cross between groups and will be referred to as links.
A single partition of the vertices a connected graph G = (V, E) into two groups V 1 and V 2 gives a graph partition H = ({G 1 , G 2 }, {l}). The link l associating G 1 and G 2 is itself a subgraph of G given by the collection of edges and the vertex set The process of repeatedly partitioning subgraphs split by this process will be referred to as recursive spectral segmentation. Empirically, we will introduce a bound on the spectral gap as a stopping condition for the recursion.

Molecular Beading
A total of three molecules, SDS, NaCl and water, are simulated in this study. The atomic to DPD bead mapping is outlined in Table S1, while each water bead consists of two water molecules. Similarly, a DPD NaCl molecule is composed of two unbound DPD beads: Na and Cl. Both of these ions are denoted to be hydrated and consist of two atomic water molecules along with the corresponding ionic atom.
The construction of the DPD SDS molecule is illustrated in Fig. S1. This linear molecule has a radial bond U (r i,j ) between each neighbouring pair of beads i and j, as per Eq. 7.
Here κ is set to 150 and r 0 is defined by the equation r 0 = (n i + n j ) − 0.01 where n i is the number of "heavy" (non-hydrogen) atoms per bead. The equilibrium bond angle U (θ ijk ) is set to 180 degrees for all angular bonds, along with a stiffness parameter κ set to 5. Figure S2: Schematic of the beading structure of SDS.

Force Field Parameters
The force-field parameters employed in this study are displayed in Table S2. S8

Mathematical Derivations of Simulation Observables
The simulation observables which we computed in the study are as follows: • The number-averaged micellar cross-sectional radius r cs .
• The number-averaged micelle backbone length l b .
steps to ensure that the number of equilibrated frames exceeded at least six of the longest observable autocorrelation timescales, and reliable statistics could be gathered. S10 A detailed explanation of this protocol is beyond the scope of this paper and may be found in the following reference. S9