Frederic Cuppens
email: frederic.cuppens@enst-bretagne.fr

Nora Cuppens
email: nora.cuppens@enst-bretagne.fr

Julien Thomas
email: julien.thomas@enst-bretagne.fr

S-TGDH, secure enhanced group management protocol in ad hoc networks

We present and analyse a secure protocol for group management, in large and dynamic ad hoc networks. The protocol we suggest relies on the TGDH protocol.

In comparison with the previous solution, our algorithm helps to uniformly dispatch the group key calculus on each node, and the global cryptographic tree is optimized. Moreover, we propose an authentication algorithm.

Our algorithm provides several well-known security properties, such as nodes authentication, messages freshness, passive attacks resistance and known key attack resistance.

I. INTRODUCTION

Group communication has been studied for a while, in wired and wireless networks. With the emergence of security needs in wireless network and more specifically in ad hoc networks (for military's applications or consumers ones, for instance with multi-players games or tele conferences), several protocols suggest a way to get secured group communications. To provide group communication privacy, we need to have a common secret key. Efficiently managing group keys in mobile networks, such as ad hoc networks, is difficult as every time the group changes (a node joins or leaves the group), the group key needs to be modified.

In this paper, we present and analyze a secured and authenticated group key management (S-TGDH), based on an existing solution, TGDH [START_REF] Kim | Simple and fault-tolerant key agreement for dynamic collaborative groups[END_REF]. The rest of this paper is organized in five sections. Section II reviews previous approaches for the distributed group key management and presents briefly the TGDH protocol. Section III and IV present our contributions. Section V gives complexity analysis of S-TGDH and some comparisons with TGDH. Finally, section VI concludes and discusses future works.

II. RELATED WORKS

A variety of solutions have been proposed for group management. Some of them rely on centralized management, for instance Simple Key Distribution Center (SKDC). The logical key hierarchy (LKH) protocol [START_REF] Wong | Secure group communications using key graphs[END_REF] is another group management protocol in which each member gets a place in a cryptographic tree managed by a centralized server which performs all the group management operations.

Other papers study decentralized algorithms. In decentralized group management protocol, each node of a group has partial keys. When we put all of them together, we are able to create the group key, needed to (de)cipher communications. Several proposals have been suggested to manage the group key. We have for example the clique suite proposed by Steiner and al [START_REF] Ateniese | New multiparty authentication services and key agreement protocols[END_REF], in which every member introduces its partial-key into the result generated by its preceding member and sends the new result to its following member.

Kim and al [START_REF] Kim | Simple and fault-tolerant key agreement for dynamic collaborative groups[END_REF] propose TGDH, which arranges keys in a tree structure (see the following section). TGDH is similar to the One-Way Function Tree (OFT) protocol [START_REF] Sherman | Key establishment in large dynamic groups using one-way function trees[END_REF] except that TGDH uses Diffie-Hellman instead of one-way functions for the group key generation.

Tree based solutions

An interesting aspect in TGDH and the other tree-based solutions is that the protocols try to dispatch the group key calculus on nodes of the group, depending on the tree structure. As our proposition relies on TGDH, we present the main aspect of this solution. More explanations are available in [START_REF] Kim | Simple and fault-tolerant key agreement for dynamic collaborative groups[END_REF].

Cryptographic principles:

The TGDH cryptographic algorithm relies on the Diffie-Hellman solution adapted to the group problem. Thus, the arithmetic operations are performed in a group of prime order p with the generator α, where p (prime integer) and α (exponentiation base) are the ones described in the Diffie-Hellman protocol [START_REF] Rescorla | Diffie-hellman key agreement method[END_REF]. As described in the following section, the notion of Group Diffie-Hellman is due to the fact that the group key is generated using the Diffie-Hellman problem.

Each member possesses two personal keys: a private partial-key K v and a public one BK v , where BK v = α Kv mod p. The public group key is a function of the other nodes public keys and the current node private key. Let's say: for v, groupkey = f (K v , BK 0 , ..., BK v-1 , BK v+1 , ...BK n ), where n is the group size and f is the group key function.

When the number of nodes increases, the number of parameters for f raises up too. However, the evaluations f (K v , BK 1 , BK 2 , ..., BK w ) and f (K w , BK 1 , BK 2 , ..., BK v ) may perform similar partial calculi (on BK 1 and BK 2 ). In order to prevent several nodes to perform similar calculi, we use the cryptographic tree to dispatch the global key evaluation.

The following example explain how the group key is computed and how the calculi are dispatched on the group's nodes. Let's suppose we have 6 members. One of trees that TGDH can produce is represented in the figure 1. One should remember that the intermediate nodes (ie the ones that are not leaves) such as <2,0> are virtual nodes, which means that they have no real existence and that their keys are computed from their children : K l,v = (BK l+1,2v+1 ) K l+1,2v mod p In this figure, the copath of the node M 1 is represented (in dark line). The copath of a member is the set of neighbour nodes (virtual or not) associated to the path we get when we want to reach the root of the tree. The important aspect is that only the BKs on the copath are needed to evaluate the group key. When we use keys of virtual nodes, it means that the BK is already a computation of children's BKs. For instance, the BK of <2,3> is a computation of those of M5 and M6. Thus, the function f (K v , ..BK 5 , BK 6 , .BK w ) is replaced by f (K v , ..BK parent , ...BK w ), which reduces the number of operations a node has to perform. By comparison, when M 1 wants to compute the group key it needs to perform f with 3 BKs, and not 5 (the other members).

In TGDH specifications, we have five membership events:

• Join: a new member is added to the group • Leave: a member is removed from the group • Merge: a subgroup is added to the group • Partition: a subgroup is insulated from the group • Key refresh: the group key is updated We describe only the two first events, due to lack of space (Merge and Partition are respectively meta-events of Join, and Leave).

Join protocol: Let's suppose we have n nodes in the group. The new member M n+1 initiates the join protocol by sending a join request message that contains its own blinded key BK. When current group members receive this message, they first determine the insertion node in the tree. If the cryptographic tree is well balanced, M n+1 joins at the root node. Otherwise, the insertion node is the shallowest rightmost node, where the join does not increase the height of the key tree. The sponsor is the rightmost leaf node in the sub-tree designed by the insertion node (see figure 2).

When the sponsor is found, it creates a new intermediate node, and promotes the new intermediate node to be the parent of both the insertion node and the new member node. The sponsor broadcasts then the new tree, which contains all blinded keys. Leave protocol: Let's suppose that a member M l leaves the group. In this case, the sponsor is the right-most leaf node of the sub-tree rooted at the leaving member's sibling node, as shown in figure 3. The former sibling of M l is promoted to replace M l 's parent node. The sponsor picks a new secret, computes all keys on its copath, and broadcasts the new set of blinded keys to the group. This information allows all members to recompute the new group key. In this section, we present the first part of our TGDH's variant: we introduce our authentication algorithm (not specified in TGDH's paper), which suggests

• a way to authenticate other nodes starting from a weak group password • a way to authenticate source message during exchanges First, we expound the main principles of the authentication part of our protocol and the corresponding authentication algorithm. Then, we present the different security properties our algorithms suite assures. The final part, about the keys management, describes which kind of key we advise to use, as it is important in security in general, though it does not immediately interact with our solution.

A. Principles

In addition to the group management problem, the authentication in ad hoc network cannot rely on a centralized authority. Thus, before creating a group, we must find a way to diffuse authentication parameters to the other nodes of the group. Often, these parameters are RSA public keys, which is the case in our proposition. However, in order to prevent from having a message number explosion, a temporary centralized authority has to be used. We call it the initiator node (M n ), as it performs the group exchange key initiation.

In many protocols, this RSA key distribution is not discussed. Thus, we can suppose that these parameters are exchanged via a private channel, or by static configuration. Theses solutions are not scalable and it is unrealistic to suppose that the static configuration can be performed in public (or personal) networks.

The Asokan-Ginzboorg protocol [START_REF] Asokan | Key-agreement in ad-hoc networks[END_REF] suggests a way to create a strong secret starting from a weak key. The figure 4 sums up the algorithm, whose secured solution is proposed by Graham Steel, Alan Bundy, and Monika Maidl [START_REF] Steel | Attacking group protocols by refuting incorrect inductive conjectures[END_REF]. In fact, it can also be used to handle RSA keys distribution: when creating a group, nodes only have to decide of a common temporary key. Starting from this weak secret key, we can use a challenge protocol (the group creation protocol) to perform an authenticated key exchange protocol. This can be seen as a strong assumption, but in fact it is quite obvious: when people want to create a group, they have a secret they all share, in order to be assured that others are the one they pretend to be (for example, they know each others' faces).

B. Algorithm

The exchange algorithm we propose is an adapted version of [START_REF] Steel | Attacking group protocols by refuting incorrect inductive conjectures[END_REF], which joins together the authentication and the group initiation processes. In the TGDH specification, no assumption is made about the group creation. It can be done for instance by the following initialization protocol: the initiator sends a create group message. Upon receiving an invitation response from a node M i , the initiator adds M i to the tree (thus, the place of M i only depends on when it sends its invitation response). After receiving all the responses, the initiator sends the group tree to each member. We have here a 3-steps protocol.

If we perform a similar protocol, by using a variant of the Asokan Ginzboorg protocol, we get the 4-steps protocol described in the table I.

I AUTHENTICATION PROTOCOL IN S-TGDH 1. Mn → M i : {publicn, c1i}p, n1, signaturen(i = 1..n -1) 2. M i → Mn : {public i , c1i, c2i} publicn , n2, signature i (i = 1..n -1) 4. Mn → M i : {c2i} public i , n3, signaturen(i = 1..n -1) 3. M i → Mn : agree, n4, signature i (i = 1..n -1)
In the first step, M n , the initiator node, asks all the nodes if they want to join the group, using the weak secret. It sends its public RSA keys and a first challenge c1i, used to authenticate the node M i . When M i receives the message, it first deciphers the message, using the same secret. It can here perform a signature check.

At step 2, M n can check if M i knew the weak secret P by testing if private n (public n (ci1 ′ )) = ci1.

The step 3 is performed by M n when it received the messages (step 2) from all the members. It can join the tree structure with this message. At this step, M i can also check M n identity, by testing if private i (public i (ci2 ′ )) = ci2.

The last step is needed for group confirmation from each node (to be sure no error occurred). If the underlying protocols are reliable, the step 4 can be ignored. But it depends on the routing protocol (for example, OLSR [CeA + 03] uses UDP). Note that each message comes with its signature, which provides identity confirmations.

We made here the assumption that the group creation is based on a common group password, as we assume that no node has several identities in the network. Thus, if this property cannot be assured, it is better to replace the common group password by a 2-nodes password.

C. Security properties

In cryptographic protocols, several security properties have to be satisfied in order to keep the confidence in the group security. We present here those properties our algorithms preserve. Most of them are taken from the clique protocols specification [START_REF] Ateniese | New multiparty authentication services and key agreement protocols[END_REF]. These different properties have to be satisfied in order to comply with two basic properties : outsider cannot read and outsider cannot send. (note : some of these properties were already guaranteed by TGDH, we mark them by T )

Group management:

• station authentication: upon receiving a message from a node, we should be able to detect if the claiming node is the real one (no masquerade attack). This is provided by the signature mechanism, using the RSA keys. • keys independence (T ): this property can be divided in two basics one, as described in the Hi-KD [START_REF] Hassan | Hi-kd: Hash-based hierarchical key distribution for group communication[END_REF] protocol :

-Forward confidentiality1 : an old node must not be able to get a group key used after it left the group -Backward confidentiality: a new node must not be able to retrieve group secret keys used before it joined the group. These two properties are ensured by the key update mechanism upon each group event and the Diffie-Hellman properties (it is impossible to get a private key, no matter which public keys an attacker has).

• non-contributory protocol: all the nodes in the group do not play the same roles, but a node's role cannot be chosen by itself. The fact that the role cannot be chosen is very important in order to prevent nodes from taking particular role. • message freshness: in order to prevent malicious nodes from performing message replays, the protocol must certify the freshness of each message. In our solution, this is assured by nonces.

Group key:

• key agreement protocol (T ): a protocol is said to be a key agreement protocol if the secret key is computed from two or more parties and that no party can predict the resulting value. to prevent the group from using the same key during a long time), prevents attackers from discovering the current session key. As all these properties are satisfied by our algorithms, we are able to keep the confidence in the group security.

D. Keys & signature management

Key validity: As we are in a distributed environment, updating a key is not immediate. Thus, it should be taken into account that an old group key should not be revoked immediately after the new one is built.

In [BBB + 06], several advices are given about key management. One aspect, the cryptoperiod, has to be taken into account for key validity. "A cryptoperiod is the time span during which a specific key is authorized for use by legitimate entities, or the keys for a given system will remain in effect." In fact, a BK (or K) key in S-TGDH has three states. The first one, pre-partial key, refers to the fact that the key as been created, but is not a part of the global secret key of each node (for instance due to the propagation delay for an update operation). Then we got the partial-key and finally the post-partial key. Some of the node use old keys (in post-partial state) and others use new keys (in pre-partial state), due to delays.

This can be summed up in figure 5. The time for each phase depends on the network, but also on the sensitivity of the information. Key size: An important aspect in cryptographic protocols is the key size. Nowadays, asymmetric key protocols use keys of 1024 bits length (and symmetric ones of 128 bits). According to a report [START_REF] Williams | A discussion of the importance of key length in symmetric and asymmetric cryptography[END_REF] of Lorraine C. Williams, 2048 bits length for asymmetric key and 128 bits length for symmetric key is far enough for the moment.

We thus recommend the use of these key lengths (depending on the country laws). For the global key, the use of a hash 128 or hask 256 functions is sufficient.

Nonces management: Two kinds of nonces can be used in cryptographic protocols. The first kind of nonces, based on random variables, assumes that all nodes know the previous used values. However, they should know all the nonces used since the RSA keys has been created, and not only since the group has been created. Thus, this solution is not well adapted to our problem. The second kind of nonces relies on temporal leashes [START_REF] Hu | Packet leashes: A defense against wormhole attacks in wireless ad hoc networks[END_REF]. This solution, most interesting for group management, implies several requirements such as :

• nodes should be able to synchronize their clocks • the propagation time can be evaluated For the first requirement, solutions as GPS can be used. Several synchronization protocols exist in ad hoc networks, but we think GPS is better due to traffic increase problems. The propagation time (or a validity time) is needed to give a validity interval for a nonce.

IV. ALGORITHMS OPTIMIZATIONS

For network protocols, two of the main complexity problems are how to manage a message and how many management messages are sent on the network. The optimizations of the first problem help to decrease the part of the node capacity associated to the protocol management, while the optimizations of the second help to reduce the consumption of energy and to increase the useful bandwidth.

TGDH proposes to create a cryptographic tree, which helps to dispatch the complexity within the nodes of the group. However, the way the sponsor is chosen and the fact that all the intermediate nodes in the tree are virtual lead to a non-uniform distribution of the global complexity, as the left children nodes would have more and more works to perform. Our idea consists in having the less virtual nodes as possible which will help to dispatch uniformly the work. The following sections describe how we update the different TGDH algorithms to remove as many virtual nodes as possible.

Moreover, having less virtual node reduces the size of the tree, which has two consequences:

• the tree propagation is faster • less BK keys are needed to build the group key In order to lessen the complexity, we present in this section some modifications of the joining and leaving algorithms of TGDH. As we modified the tree aspect (by removing as many virtual nodes as possible), several cases have to be studied when these events occur. In order to satisfy the station authentication property, each message is sent together with its signature, using the RSA keys discussed before.

A. Joining algorithm

Join authorization: As we have now a secured and fully authenticated group, adding a node to the group has to be authorized. We make the assumption that when a node that belongs to the group wants another one to join the group, the authorization is valid. As the group key is known by all the members, this assumption is not so strong: nothing can prevent a node from giving the group key to a node which does not belong to the group. However, this assumption does not take into account malicious behaviors for which the goal is to attack the group management protocol, and not only to listen to communications. This problem is discussed in the paragraph Trustworthy join, as we can consider, in our case, multiple malicious nodes as multiple identities of a single malicious node.

Using the terminology of TGDH, the node that invites the new node is the sponsor. It has to perform three actions: first, it updates its partial keys, in order to assure the backward secrecy property. Then, it exchanges information with the new node. Finally, it sends the join message to all the members of the group, acting as the new node's sponsor.

Algorithm principles: When a sponsor adds the node to the tree, it should do it in a way that reduces the number of virtual nodes. One choice would be to get the best place in the whole tree, and to assign this position to the new node. However, this is not possible in a large network, or in an active ad hoc network, as conflicts could occurred (several sponsors would choose the same position, according to the fact that message propagation is not instantaneous). We then should perform local calculi, as there would not have any conflict.

Depending on the local context, two different joins can occur:

• if the parent is a virtual node, then the current node can (see the following requirements) choose to replace it by the new node (figure 6). The key replacement algorithm is then simple: the BK key of the new node is computed by the sponsor node and propagated with the new subtree in the ADD message. • if the parent node is not virtual, the node create a subtree with a virtual parent, whose children are itself and the new node (figure 7). Then, the sponsor node creates the virtual node BK and propagates it. The key update algorithm (which comes after the BK keys replacement) is the same as the group creation one.

Algorithm requirements: In order to optimize the group management process, the algorithm must meat the following requirements:

• only leafs can act as sponsor, in order to avoid merges • only left child can replace parent virtual node, again in order to avoid merges • the global tree must stay well balanced, to dispatch work uniformly. This means that a node can create a sub-tree only if its depth does not exceed too much the average one. These requirements are quite strong, but not so difficult to implement: when a node which is not a sponsor wants to add a new node, it can delegate it to another node. The fact that all the messages are authenticated helps to do that. Trustworthy join: In order to have correct group decisions, it is important to prevent a node from having multiple identities. If this is not possible, we propose a modification of the S-TGDH joining algorithm: a node is allowed to join the group only if at least τ A nodes (similar to the threshold cryptography [START_REF] Stanis | Adaptively secure threshold cryptography: Introducing concurrency, removing erasures[END_REF]) allow it to join the group. The authorization to join a group can depend on several parameters, such as reputation [START_REF] Liu | Enhanced reputation mechanism for mobile ad hoc networks[END_REF].

Thus, the new node has to collect proofs from τ A nodes, where proof = {nonce, nonce private RSA key }.

This threshold value τ A depends on several context properties such as the trustworthiness of each node and the dispersion of the nodes (scalable solution).

B. Leaving algorithm

In TGDH, a node departure is not authenticated. Thus, a malicious node can eject people from a group by only sending fake message, taking the identity of the node it wants to eject. With S-TGDH, the message authentication guarantees that this is impossible. Thus, the leaving algorithm is similar to the TGDH one. But we also introduce a leaving proof, to improve the security. Finally, we discuss about the silent leaving problem.

Leaving sponsoring: In order to properly leave a group, a node must tell its sponsor node it wants to leave. The sponsor then updates its group keys in order to respect the forward secrecy property. However, the way the sponsor node is chosen is not trivial. Depending on the way we do it, we could obtain conflict between several leaving nodes, but also with new nodes (which is more complicated to deal with).

According to how we manage the sponsor node for a new node, the election of the leaving sponsor node has to happen in the following way (see figures 8 & 9 ): the sponsor node is the rightmost leaf in the neighbour subtree. The use of its neighbour and not its own children prevents merge with local adding (as adding a new node only modifies the subtree of the sponsor node). The use of the right child provides a deterministic way to choose a sponsor, to perform check operation (it could have been the left, but using the right one dispatch the work well).

As for the adding operation, the leaving operation generates a BK message from the sponsor, which in turn generates a global key modification.

The leaving of a node generates a hole in the tree. In order to prevent from having an increase of virtual nodes, several cases should be taken into account:

• the parent node is virtual (figure 8). Then, the neighbour node can go one step upper in the tree, taking the place of this virtual parent • the parent node is not virtual (figure 9). Then, a virtual node takes the place of the parent node, and this parent node takes the place of the leaving node Our algorithm provides an efficient way to realize leaving operation without merges. However, we can see a conflict problem when a node and its sponsor want to leave at the same time. In order to prevent deadlocks from happening, due to partial view problems, nodes have to wait for their sponsor's confirmation, before leaving the group. Leaving proof: Though the messages are authenticated, fake leaving can happen from the inside : the sponsor node can claim to have received a leaving message from a node, and then every node in the group can trust it. Thus, when a node begins to act maliciously, it can destroy the group easily.

In order to make the group more stable, we decide to implement a proof system : when a node wants to leave the group, it sends a leaving message to its sponsor, message in which it puts a proof. The sponsor can then check this proof and forward it to all the members of the group. This proof uses the signature system {p, (p) private RSA key }.

Silent leaving: As we have seen, deadlocks could have occurred if two nodes want to leave at the same time. However, other deadlocks can happen if nodes do not comply with the protocols. This can be intentional, or not (for example due to crash problems).

If this happens, local nodes can detect it and inform the others that a bad behaviour occurred. However, if only a single node is needed to correct this problem, a malicious node can here again send fake messages and destroy the group. Resolving this kind of attacks is quite complex, as the malicious node already belongs to the group. A convenient way to have a scalable and trustworthy solution is to use threshold system, in a similar way as threshold cryptography [aJL00]: a decision is taken if τ L nodes agree on it. The value of the threshold τ L depends on several parameters, such as the ones stated for τ A . Obviously, the bad behaviour detection is efficient if less than

• τ L nodes are contaminated, to have correct alerts • n -τ L nodes are contaminated, to signal bad behaviours Using reputation mechanism can also prevent malicious nodes from performing collusions: only nodes with a good reputation are allowed to participate in the group decision process. By putting together the leaving and join restrictions, we get a more resilient group management protocol.

V. COMPLEXITY ANALYSIS

A. Optimizations analysis

The figure 10 sums up how a node can manage BK messages in the TGDH protocol, which are the most important messages exchanged when the group is created.

One can clearly see that the protocol is simple, but the problem with TGDH is that a lot of BK messages for virtual nodes are sent. The number of messages sent when a node updates its BK depends on its depth in the tree. Let's consider n nodes in the group. In TGDH, the depth of the tree is given by the formula h = E(log2(n))+1 and h BK keys have to be rebuilt (in order to get the tree root). Another problem is that as all intermediate nodes are virtual, so one of the leaves has to propagate these messages. The average number of BK keys a node should manage is

hs-1 i=1 (i • (1/2) i+1 ) + h s • (1/2) h
. The optimization we propose with S-TGDH (illustrated by the figure 11) helps to reduce the number of messages sent during a BK key update process, and also the number of messages sent by each node. As for TGDH, the number of messages sent during a key removal depends on the depth of the considered node. But contrary to TGDH, all the intermediate nodes are not obligatory virtual. Ideally, we have no virtual node. Thus, the average number of sent messages is h s = E(log2( n 2 + 1)) and the other interesting part is that a node may have to take into account only one BK message, if the parent is not a virtual node. We also notice that the number of messages associated to a BK removal process depends on an implementation choice: either we can use a multicast process if the routing protocol allows it, or we can use unicast messages. Depending on this choice, the number of messages sent by the source and the intermediate nodes is quite variable. Thus, as it does not change the comparison between TGDH and S-TGDH, we choose to compare the number of BK keys that have to be changed and the average per node.

B. Algorithms complexity

In order to estimate the algorithm cost of each subprotocol, we make the assumption that an optimal messages propagation algorithm can be chosen (see previous section).

The results are summed up in the tables II and III.

Group creation:

The number of messages sent during the creation process (ie the authentication phase) is 4 • (n -1), where n is the number of nodes if we consider unreliable network communications (3 • (n -1) otherwise). Comparing to the 3 • (n -1) messages suggested by TGDH, the process produces an increase of 33% in the worst case though our approach provides a fully authenticated protocol. Besides, as explained below, this overhead is produce only during the group creation process and is compensated by the optimisations made in the other processes. All the algorithm complexity is supported by the initiator (as other nodes only have to check the challenge from the initiator node and to decrypt messages). For the initiator node, the complexity is linked to the problem where to insert a node in the graph. With the implementation we propose, the complexity is minimal: the maximal cost is obtained if no virtual node is caught during the sub-tree analysis. We then have h • single cost. The global cost is n i=0 E(log2(i+1)) which give a complexity of n•log2(n).

Group initiation:

In this case, the number of nodes that get a BK message is maximal, since all the leaves send their BK keys. So we have n leaf = 2 Hs-1 = n+1 The different algorithm costs are evaluated according to our implementation of S-TGDH. Moreover, the algorithm complexity in S-TGDH is quite the same as in TGDH as the additional checks on S-TGDH only happened on the copath, which makes them fast. Updating a key: The number of nodes that get a BK message is Hs i=Hs-hs+2 (2 i -2) where H s = E(log2(n + 1)), as mentioned before (the sum on i is due to the removal of the intermediate BK keys, until we reach the tree root). The cost for the K messages is h s • Kcost.

C. Authentication overhead

With S-TGDH, all the messages are authenticated. The overhead is produced by the signature attached to each message, which is in fact given by the formula signature = hash((message) key ). The hash function, as for other cryptographic protocols, generates an output of constant length. For instance, the commonly used md5 hash function produces an output of 128 bits.

By default, the packets maximum size is the same as for all the Ethernet protocols: 1500 bytes. If we consider a transfer that tries to use as many payloads as possible, we got a decrease of 3% of the performance (using a 256 bits length signature). But very often, the packets do not have the maximum size. Thus, the performance decreases come only from the time spent to send the packets. For the WiFi, it takes 250 µs with 1 mbps (slowest speed) to send the signature.

Another impact of the packet's encryption is the time needed to decrypt it. Fortunately, only the recipient has to decrypt the packet, as the addresses are not encrypted. As we use a symmetric cryptographic process, the deciphering process is quite fast (for instance, one can have to perform only a xor operation).

D. Graphical comparisons

The following graphics illustrate the different results we described in the previous sections and compare them with the one we got for TGDH. We can see that both protocols produce the same cost (with better performance for S-TGDH), but S-TGDH generates fewer messages and distributes equitably the cost on each node.

VI. FUTURE WORKS

We have studied how to enforce security in the TGDH protocol, to propose S-TGDH, which improves several security aspects in group management. Moreover, we propose an authentication algorithm embedded in the S-TGDH group creation protocol. The optimization we propose helps to dispatch uniformly the complexity over every node. However, message confidentiality in group communication is not the only security aspect to study. Thus, in future works, we will study other aspects such as anonymous routing and group cooperation, which would be interesting extensions of our current solutions. Moreover, proposing a way to evaluate reputation at the group layer (in addition to the ones at the lower levels) is useful to have a good group decision, in order to detect internal attacks (see section IV-B). Furthermore, a lot of security flaws in security protocols come from the fact that these protocols were not formally specified. So, we are currently investigating formal modeling and proofs of our current proposition. Finally, the current implementation of S-TGDH is a java simulator, which simulates the algorithm behaviour. An important dimension in ad hoc network is the mobility [START_REF] Davies | Evaluating mobility models within an ad hoc network[END_REF]. Thus, implementing our protocol using the well known simulator NS-2 will help us to test our solution in mobile networks, using adequate mobility models, and to analyze the interaction with the underlying layers (for instance, OLSR or AODV).

Fig. 1 .

 1 Fig. 1. TGDH tree example

Fig. 2 .

 2 Fig. 2. TGDH -join protocol

Fig. 3 .

 3 Fig. 3. TGDH -leave protocol

Fig. 4 .

 4 Fig. 4. Improved Version of the Asokan-Ginzboorg protocol

  Fig. 5. Cryptoperiod explanation

Fig. 7 .

 7 Fig. 6. Adding a node: first case

Fig. 8 .

 8 Fig. 8. leaving operation: first case

Fig. 9 .

 9 Fig. 9. leaving operation: second case

  Fig. 10. TGDH management messages

Fig. 11

 11 Fig. 11. S-TGDH management messages

  H s = E(log2(n + 1)). The global cost is then n leaf • Hs i=2 ((2 i -2) • BKcost) for the BK messages cost, and (2 hs-1 -1) • Kcost for the K messages cost where h s = E(log2( n 2 + 1)) is the mean depth for the node to take into account. (for TGDH, only the leaf receives BK messages , so cost = n • Ht-1 i=1 ((2i -1) • BKcost))

•

  Kcost: SC • O(log2(n/2)), mainly in order to build its public key • BKcost: SC • O(n) + SC • O(log2(n/2)), mainly to create the parent private key, and perform checks SC (Security Cost) represents the cost to create the key.

TABLE

  • implicit key authentication(T ): each member M i of the group M is assured that no member M q / ∈ M can learn the group key unless it is disclosed by a dishonest member M j ∈ M • key integrity (T ): the group key is computed only by using values issued by indoor members of the group. Extraneous contribution(s) to the group key cannot be tolerated even if it does not provide the attacker(s) with any additional knowledge.

• Perfect Forward Secrecy (PFS) (T ): the compromise of a key (long-term or session's key) cannot result in the compromise of past session keys.

• known key attack resistance: compromise of session keys does not allow: 1) a passive adversary to compromise keys of other sessions, or 2) an active adversary to impersonate one of the protocol parties. The signature mechanism prevents an adversary to impersonate a node, unless it manages to obtain this node's RSA private key. • passive attacks resistance: a protocol is said to be resistant to passive attacks if it can prevent people from discovering secrets (such as the group key) by passive attacks, like eavesdropping. The group encryption mechanism, associated to group key updates (in order

As precised before, this definition of the Forward confidentiality is a Hi-KD definition[START_REF] Hassan | Hi-kd: Hash-based hierarchical key distribution for group communication[END_REF], and not the international definition, which is equivalent to the PFS definition mentioned later. We add this Hi-KD definition as we think that it defines well the security aspect of the group update process.


Maria Panoukidou

Simon Weir

Yair Valerio Sorichetti

Gutierrez Fosado

Martin Lenz

Davide Michieletto
email: davide.michieletto@ed.ac.uk

Yair Gutierrez Fosado

Runaway Transition in Irreversible Polymer Condensation with Cyclisation

published or not. The documents may come   L'archive ouverte pluridisciplinaire

The process of polymer condensation, i.e. the formation of bonds between reactive end-groups, is ubiquitous in both industry and biology. Here we study generic systems undergoing polymer condensation in competition with cyclisation. Using a generalised Smoluchowski theory, molecular dynamics simulations and experiments using DNA and T4 ligase, we find that this system displays a transition, from a regime with finite-length chains at infinite time and dominated by rings to one dominated by linear polymers that grow in time. Finally, we show that fluids prepared close to the transition may have profoundly different compositions and rheology at large condensation times.

Linear polymer condensation is the process by which two polymeric end groups react and form a bond. It is one of the most common industrial processes to create long linear polymers [1], and is also relevant to DNA repair [2]. Classic theory of linear polymer condensation assumes that, in absence of loop formation, the average length of the reacting polymers diverges as l = 1/(1-p) where p is the extent of reaction [1,3]. However, loop formation is expected to be favourable in certain conditions [4][5][6][7], yet difficult to precisely control during the synthesis [8]. The formation of loops in polymer condensation poses open challenges to both the theoretical and experimental polymer communities and significant efforts were made to understand the ring-chain equilibrium in reversible systems [7][8][9][10][11][12]. During irreversible polymer condensation, the formation of loops is expected to have a two-fold impact: structurally, it reduces the amount of reactable substrate that can form longer chains and, physically, it affects the rheological properties of the fluid [13]. Indeed, ring polymers are among the least understood classes of polymer architecture and have nontrivial static [14,15] and rheological signatures [13,[16][17][START_REF] Michieletto | Proc. Natl. Acad. Sci[END_REF][START_REF] O'connor | [END_REF][20]. Due to the importance of chain topology in both dilute [21] and dense [13,[START_REF] Michieletto | Proc. Natl. Acad. Sci[END_REF][START_REF] O'connor | [END_REF][22][23][24][25][26][27] polymer solutions, it is thus important to understand the assembly kinetics of ring and linear polymer topologies during polymer condensation.

Here we study irreversible linear polymer condensation using a combination of theory, simulations, and experiments. First, we find that the competition of polymer condensation and cyclisation can be well described by a modified form of the Smoluchowski coagulation equation [28] with a sink term that captures the formation of loops. By spanning a broad range of polymer concentrations c, we discover that there is a runaway transition characterised by the fact that linear chains grow quickly enough to permanently escape cyclisation. This transition separates a regime in which all the polymers are converted into rings at infinite time (and hence do not grow with the extent of reaction), from one in which the mass fraction of ring polymers is finite, and the length of the linear chains diverges in time. Importantly, and in contrast to typical textbook pictures, this transition occurs at a finite, non-zero looping probability.

Finally, we show that as a consequence of this transition, systems prepared at similar initial concentrations close to the runaway transition and undergoing irreversible condensation will display markedly different architectural and rheological features at large enough condensation times. We conclude our paper by discussing the implications of our findings in the design of soft materials and molecular cloning.

Results -We first simulate linear condensation by means of Molecular Dynamics. We use a bead-spring polymer model [29] to simulate l 0 = 174 beads-long chains (bead size σ) diffusing in a box with periodic boundary conditions at monomer concentration c (see SI for details). The condensation process is modelled by creating an irreversible bond between two end-beads that are closer than 1.1σ in 3D space. The condensation is stochastic and occurs every 10 3 τ B (τ B = γσ 2 /k B T is the Brownian time and γ the friction) with probability p L = 0.1 (Fig. 1a,b). During the simulation, we take snapshots of the system and record the list of bonds to reconstruct length and topology of the polymers (Fig. 1c). Over the simulation time, the number of initial linear chains decreases due to the formation of (i) longer linear polymers or (ii) circular chains (Fig. 1d). Additionally, lower monomer concentrations c promote the formation of more rings at large times and a slower decrease of the linear species. We also note that (i) the number frac- tion of rings converges to a finite value at large time, and that (ii) while the number of linear chains appears to go to zero, their length increases (as shown in detail below). Accordingly, the (number) average length of polymers, i.e. l = l 0 i n i i/ i n i , grows more quickly for larger c (Fig. 1e). Thus, we conclude that loop formation competes with the growth of the chains, and that cyclisation is dominant in dilute systems.

To validate our simulations we perform irreversible condensation on linear DNA using T4 ligase (NEB). This is an enzyme that consumes ATP to form a covalent bond between two proximal and complementary doublestranded DNA ends. More specifically, we perform irreversible condensation (also referred to as ligation henceforth) on a monodisperse solution of linear, l 0 = 6.5 kblong, DNA at low concentrations c = 2 ng/µl 0.01c * (c * = 3l 0 M w /(4N A πR 3 g ) 0.2 µg/µl is the overlap concentration, M w = 650 g/mol the molecular weight of a DNA basepair and N A the Avogadro number). After adding T4 ligase, we draw aliquots at time intervals and heat-inactivate the reaction by heating the aliquot at 65 • C for 5 minutes. We then split the aliquot and treat one of the two sub-aliquots using exonuclease (RecBCD, Lucigen), an enzyme that digests linear, but not circular, DNA. Finally, we treat all aliquots with nickase, to relax torsional stress and remove supercoiling from the looped DNA [START_REF] Bates | DNA topology[END_REF]. By running the aliquots on agarose gels we can visualise and compute the fraction of molecules in the linear and ring, monomeric, dimeric, etc. states. Fig. 1f reports a picture of one such gel, showing a single band of monomeric linear DNA (as it disappears after exonuclease treatment) at t = 0, evolving into three bands, one of which is exonuclease resistant (a monomer ring) at larger times. In Fig. 1g we plot the relative abundance of these populations, from which we obtain the number average molecular length l(t) (Fig. 1h), which plateaus towards a finite value at large ligation times.

To understand these measurements we propose a modified Smoluchowski coagulation equation:

dn l dt = 1 2 i+j=l k 1 (i, j)n i n j - i k 1 (i, l)n i n l -k 0 (l)n l (1)
where n l is the number density of chains of length l. The first two terms describe polymer condensation while the latter captures cyclisation. The rate of condensation k 1 is assumed to take the form [START_REF] De Gennes | [END_REF]32]

k 1 (i, j) = κ 1 i -α + j -α (i ν + j ν ) ∼ (D i + D j )(R i + R j ) (2
) where D i and R i are the diffusion coefficients and size of polymer of length il 0 , with α the dynamics exponent describing the dynamics and ν the metric exponent capturing the length-dependence of their size. The cyclisation rate is taken to be k 0 (l) = κ 0 l µ , where µ = -4ν. Note that this is different from the classic Shimada-Yamakawa theory [33,34] which would predict µ = -3ν at l l p (l p the persistence length) because we (i) are out-of-equilibrium and (ii) account for the subdiffusion of the polymer end within the volume of the coil (see SI). Since we initialise our simulations and experiments below entanglement conditions, we fix α = 1 as expected for Rouse dynamics and ν = 0.588 as expected for selfavoiding polymers [START_REF] Doi | The theory of polymer dynamics[END_REF]. We numerically solve Eq. ( 1) and fit the data of mean length versus time, l(t) , obtained in simulations and experiments via the free parameters κ 1 and κ 0 (Figs. 1e,h).

The key parameter in our system is the ratio of the rates at which linear polymers and rings are generated. We thus define a dimensionless "cyclisation parameter" κ ≡ 2κ 0 /(n 0 κ 1 ), where n 0 is the number density of monomeric chains of length l 0 at the start of the simulation or experiment. Albeit related to the classic j-factor In the network representations each circle i represents a monomer and its colour represents the number of chains connected to i: 0 (grey), 1 (cyan), or 2 (blue). A linear monomer is grey, a ring monomer blue, a linear dimer has two cyan nodes, and a dimer ring would have two blue nodes. e Fraction of monomers in the giant connected component as a function of κ. f. Average polymer length (including both linear and ring) at a large time as a function of κ. g. Average length of only linear chains at a large time as a function of κ. h. Average total number of monomers belonging to rings Mrings = Nrings × lrings(t 1) at large simulation time. Note that some some simulations may not display any rings, in turn bringing the average Mrings below 1. MC = "Monte Carlo"; MD="Molecular Dynamics".

employed in DNA looping [5,33], our cyclisation parameter is more naturally interpreted as the number of rings formed for every two linear chains that are fused together. Intuitively, this number determines the final topological composition of the system. At κ 1, we expect the final state of the system to be dominated by rings, while for κ 1 to be dominated by linear chains. Importantly, since k 0 ∼ l(t) -4ν the probability of ring formation decreases in time as the average length of the linear chains increases. Accordingly, and in spite of the fact that our system has a ring-only irreversible absorbing state, we conjecture that the strongly decreasing looping probability may effectively yield a very long time-transient in which the system is dominated by entangled linear chains with circular contaminants.

Since we expect the Smoluchovski equation to be valid only in the limit of low density where three-body interactions are negligible, the values of κ 0 and κ 1 should be independent on concentration only when c is small enough. By plotting κ ≡ 2κ 0 /(n 0 κ 1 ) as a function of c/c * (where c * is computed at the beginning of the simulation or experiment) we show that κ scales as n -1 0 ∼ (c/c * ) -1 in both simulations and experiments until c c * where it starts to deviate (Fig. 2a); this confirms that the Smoluchovski approximation is valid in this range of concentrations. Importantly, in Fig. 2a we also identify the crossover value κ = 1 (at which the initial cyclisation rate is larger than the dimerisation rate) around c/c * 0.1 -0.2. We note that the agreement be-tween simulations and experiments is excellent for small c/c * . Quantitative analysis of gel electrophoresis images at larger c/c * is challenging due to the poor separation of multimeric bands.

Runaway Transition -The results in Fig. 1 suggest that at large c/c * the chains tend to grow longer, thus suppressing cyclisation; at the same time, the density of reactive ends and the speed of spatial exploration of the chains become smaller, thus suppressing dimerisation. Due to this kinetic competition, we ask whether the system can truly display a "runaway" phase, defined as a regime where at least one chain permanently escapes cyclisation and diverges in time. In Fig. 2b-c one can visually appreciate that at large reaction time, rings (blue) are abundant at low c/c * while linear chains (grey) are more abundant at large c/c * . Using a network representation (where monomer are converted into nodes and ligation events into edges) we also show that these systems display qualitative different network topologies (insets of Figs. 2b-c). At small c/c * (large κ) the network of monomers is mostly disconnected; accordingly, even when the fraction of unreacted bonds goes to 0 at t → ∞, the average length of the polymers does not diverge. On the contrary, at larger c/c * we observe only few rings and some very large chains with very low cyclisation probability. A natural quantification is obtained by measuring the fraction of nodes in the giant connected component (GCC), i.e. the largest cluster of connected monomers corresponding to the longest chain in the system. At Since MD simulations are computationally expensive, we employ Direct Simulation Monte Carlo (DSMC) [START_REF] Garcia | [END_REF][37][38] to solve the Smoluchowski equation and reach larger system sizes up to 10 5 chains (see SI). In Fig. 2f, we plot the number averaged chain length at an arbitrarily large time when the DSMC code has evolved the system as long as possible and has generated only a single linear chain accompanied by a number of rings. This figure suggests that the linear-dominated regime (κ < 1) displays an average polymer length at large reaction time that scales as l(t 1) ∼ κ -1 ∼ c/c * . At the same time, the number average length of the linear population appears to reach a plateau that scales with the system size (Fig. 2g), indicating a runaway transition at κ 1. Additionally the fraction of mass "lost" in forming rings decays as f ring ∼ κ and is thus negligible for small enough κ (Fig. 2h).

Although our simulations support the notion that a low value of κ will result in linear chains of increasing lengths and vanishing cyclisation rate, they are fundamentally limited to finite-size systems where the cyclisation rate of the largest chain never rigorously goes to 0. To firmly establish the existence of our runaway transition, in the SI we use an asymptotic theory to analytically demonstrate that the mass fraction of linear polymers goes to a finite limit at t → ∞ in a thermodynamic system. We find that the key condition to ensure the runaway transition to exist is that the cyclisation rate k 0 = κ 0 l µ decays strongly enough. More specifically, we require the exponent µ to be µ = -4ν < -4α/5 or ν > α/5. This condition is always met in the Rouse unentangled (α = 1) regime, and also in the entangled (α = 2) regime, provided that the polymers are not fully collapsed (ν = 1/3).

Dynamics -To test the consequences of the runaway transition on the dynamics, we perform microrheology experiments (see SI for details) where we track the dynamics of 800 nm PVP-coated polystyrene beads added in a solution of DNA that has been treated with either 40U T4 ligase for a week (and thus to full extent of reaction) or with buffer for a week (control). At low concentrations, we observe that the mean squared displacement of the tracer particles M SD(τ ) = r(t + τ )r(t) is unaffected by the ligation of the DNA (Fig. 3a). On the contrary, for c/c * ≥ 0.1, we see that the MSDs of the tracers in the ligated systems are much slower and displays a stronger subdiffusive behaviour than the control (Fig. 3a). From the MSD, we extract the large-time diffusion coefficient D of the tracers and in turn the effective viscosity of the sample via the Stokes-Einstein equation [START_REF] Hansen | Theory of simple liquids: with applications to soft matter[END_REF]. The plot of the normalised viscosity (Fig. 3b) suggests that a dynamical transition takes place around c/c * = 0.1 -0.2 (or κ 1 -2) which matches the structural runaway transition seen before (Fig. 3e-h).

After the transition, the viscosity increases exponentially with the concentration (see inset of Fig. 3b). This suggests a relaxation process dominated by endretraction [START_REF] Doi | The theory of polymer dynamics[END_REF], possibly due to the threading of very long linear chains through small rings [13,20,[START_REF] Roovers | [END_REF][41][42] or through pseudo-knotted parts of their own extremely long contour [43,44]. We note that, especially at large c/c * , the ligated solutions extremely elastic and the passive tracers do not display a freely diffusive behaviour even after a lag time of ten minutes. We thus argue that the reported η/η 0 may be lower bounds at large c/c * , which would render the transition even more dramatic. All this implies that, intriguingly, near the transition c/c * 0.1, systems prepared at similar concentrations may display extremely different rheology at large condensation times.

Conclusions -We have studied a system of chains undergoing irreversible condensation in competition with cyclisation. We have shown that the key adimensional parameter controlling growth kinetics is κ = 2κ 0 /(n 0 κ 1 ); naturally interpreted as the number of rings formed for any one dimerisation. At large concentrations (or κ < 1) dimerisation is kinetically favoured and drives the growth of linear chains. While growth disfavours cyclisation, it also reduces the number of available reactive ends and the annealing rate of the chains (see Eq. ( 2)), disfavouring further growth. In spite of this, we discover that the net result of this kinetic competition is a runaway transition for κ < 1 if the cyclisation rate decays strongly enough with polymer length, i.e. with ν > α/5. In these conditions, the fraction of monomers transformed into rings is finite and small, thus leaving the rest of the monomers available to form a permanently growing linear chain. Finally, we discover that the runaway transition has deep consequences on the rheology, and triggers an exponential increase for κ < 1 (or c/c * > 0.1).

Our results suggest that it may be possible to tune the final topological composition of ligated systems by judiciously choosing c/c * . For instance, the most likely regime to form large rings and ring-linear blends [13,41] is near the transition c/c * 0.1. Mixing polymer families with different reactive ends further enhances the designability as it introduces different c * for each family. Interestingly, our results also suggest that Olympic gels [22,45,46] are unlikely to form during homopolymer condensation as they would require c/c * 1 where loop formation is not favourable. Our results can also inform more efficient ways to perform genomic engineering, e.g., cloning of plasmids ought to be done at c/c * < 0.1 whereas synthetic chromosome assembly at large c/c * . Finally, it may be possible to couple dissipative DNA breakage reactions [47][48][49] with ATP-consuming ligation to create self-sustained and non-equilibrium solutions of topologically active polymers.

METHODS

Molecular Dynamics Simulations

Force-Field and equilibration process

We model a 6,500 bp-long linear DNA molecule as a bead-spring polymer made of l 0 = 174 beads. The total number of polymer chains is N m = 200. The polymers are modelled via the Kremer-Grest model [1]. Each bead has a diameter σ = 13 nm (or ∼ 38 bp), modelled as a truncated and shifted Lennard-Jones potential (WCA)

U LJ (r) = 4 (σ/r) 12 -(σ/r) 6 + 1/4 , (1)
for r < r c = 2 1/6 σ and 0 otherwise. Here r represents the distance between beads and = 1.0 parametrises the strength of the potential. The radius of the bead, σ, defines the length units in our system. Consecutive beads are connected through a permanent Finite Extensible Non-linear Elastic (FENE) bond

U FENE (r) = -0.5KR 2 0 log 1 -(r/R 0 ) 2 (2)
with K = 30 /σ 2 and R 0 = 1.5σ. The bending stiffness of the polymer is controlled by a Kratky-Porod interaction

U b (r) = k B T l p σ (1 -cos θ), (3)
which regulates the angle (θ) defined by the two tangent vectors connecting three consecutive beads along the polymer. Here, l p = 4σ = 150 bp is the persistence length of DNA. The solvent is simulated implicitly using a Langevin thermostat, so that the time evolution of our system is governed by the stochastic partial differential equations mr = -ξ ṙ -∇U + 2k B T ξη (4) * joint first author † corresponding author, davide.michieletto@ed.ac.uk where r is the position of a particle, ξ its friction, m its mass, U the sum of the interaction potentials discussed above and η white noise with unit variance. The diffusion timescale is τ B = ξσ 2 b /k B T . The integration of the Langevin equation is done with a velocity-Verlet algorithm, using a time step ∆t = 0.012τ B in LAMMPS [2].

Various monomer densities were considered, ranging from 10 -2 c * to 1c * , where c * = 0.012σ -3 is the monomer concentration at which the polymers start to overlap. The overlap concentration c * was measured by computing the radius of gyration R g of the polymers in equilibrium at very large dilution (c * = l 0 /(4/3πR 3 g )). All the systems were equilibrated for a sufficient amount of time to ensure that the polymer chains have moved at least a distance equal to R g .

Topology reconstruction

After the equilibration step, 40 replicas of production runs were started for each monomer number density considered. The ligation is performed stochastically, and is attempted every 100 time steps between two polymer end beads that are closer than R c = 2 1/6 σ using the fix bond/create LAMMPS command. The probability of successful ligation (i.e., bond formation) is set to f = 0.1. Once ligated, the bond formed between the polymers is irreversible and cannot be broken, therefore accounting for the formation of a covalent bond between the DNA fragments. During the ligation process snapshots of the system are taken every 10 6 time steps on both the 3D coordinates of the beads and the bond list at those time steps. From the bond list we can, later on, reconstruct the topology of the individual polymers, i.e. if fused with others to form linear chains or if circularised.

Once the production runs finished, the trajectories were analysed and the topology was reconstructed using our python code (https://git.ecdf.ed.ac.uk/ taplab/dna-ligation.git). We will refer to this code from now on as TR code. The code takes as input the instantaneous trajectory and bond list from LAMMPS and checks for newly formed linear and ring chains. The output of the python code is a file containing the number and length of linear chains that have formed in a arXiv:2210.14010v1 [cond-mat.soft] 25 Oct 2022

given simulation time step. Similar files are produced for the ring chains. These files are then used to calculate the average length and the number of linear/ring chains figures.

The starting point of the topology reconstruction algorithm is an array b of size N b ×2; each row b i = (id 1 , id 2 ) represents the IDs of atoms that are bonded within the system:

b =      id 1 id 2 id 3 id 4 . . . . . . id Nn-1 id Nn      ,
Since not all particles are linked together, some do not appear in the array b. In order to avoid operations with large sparse arrays, the matrix b is mapped to b that contains only indexes from 1 to the maximum number of atoms connected M : b → b.

The next step of the TR algorithm is to create a connectivity matrix C based on the list b. Each row of C represents an atom index and consists of three components C(id i , :) = (id i-1 , id i+1 , f lag). Since in our case a particle can be linked with two more particles the first two components of each row id i-1 , id i+1 represent the connections of particle id i (note that id i-1 and id + i + 1 are not necessarily consecutive in 1D but can be any other particle bonded to particle i). The third component, f lag, takes only the values 0, 1 and accounts for the particles id i that already belong to a polymer chain. The flag column of C is initialised to zeros. During the reading process of the connectivity matrix the algorithm switches the flags to 1 of the particles that are already considered to belong in a chain. Rings are extracted in the same manner and a ring is found if the current atom index is the same as the starting atom index. This reading process outputs N c arrays which have different lengths and each of them contains the particle (mapped) ids that are connected in a polymer chain. The final step of the TR algorithm is to map the atom indexed back to the original ones: M -1 : b → b. This algorithm is very generic and can be applied also in cases where the atoms are not initially in polymers as in our case, but rather individual atoms that can connect with each other during the simulation.

Experiments

We consider a 6500 bp-long plasmid (referred to as 1288 here) which is converted into linear form by using a restriction enzyme (XhoI). This linearisation step is checked on a gel electrophoresis. The equilibrium radius of gyration of this linear DNA molecule is about R g l p L/3l p 0.2 µm (in agreement with diffusion data from Ref. [3]). This yields an overlap concentration c * = 3M/(4N A πR 3 g ) 0.2 mg/ml where N A = 6.023 10 23 1/mol is the Avogadro number and M = L × 650 g/mol is the molecular weight of the plasmid. For the low DNA concentration experiments we set ourselves at 0.01c * , i.e. c = 2 µg/ml. To perform ligation we use T4 ligase (NEB, M0202L, 1U corresponds to 0.5 ng or 0.00735 pmoles of protein according to Ref. [4]), and work at 1x T4 ligase reaction buffer concentration, which contains 1 mM ATP.

To classify the topology of the DNA under ligation, we perform time-resolved gel electrophoresis. We prepare a master solution of DNA at desired concentration, 1x ligase buffer and 2 U/µl T4 ligase. Aliquots are taken at different elapsed times, the T4 ligase is heat-inactivated at 65 • for 15 minutes. Aliquots are treated with Nb.BbvCI Nickase (NEB, R0631L) to relax the supercoiled population. A subaliquot is created and treated with exonuclease (PlasdmiSafe, Lucigen) to remove the linear population. The resulting aliquots are run on a gel: we load 20ng of DNA from each aliquot onto a 1% agarose gel prepared using 1x TAE buffer. A standard λDNA -HindIII digest (NEB,N3012S) marker is also loaded. The gel is run at ∼ 2.5V/cm for 5 hours and post-stained with Sybr-Gold (ThermoFisher) for 30 minutes. A Syngene G-box and Genesys software is used to image the gels.

The combination of nickase (relaxing the DNA supercoiling) and exonuclease (fully digesting linear DNA molecules) allowed the topology of the DNA in each band to be unambiguously identified. Further, the λDNA -HindIII digest marker confirmed the bands were of the correct size for monomer and dimer lengths. To extract the relative amount of molecules in each lane we compute, using ImageJ, the intensity of each lane and account for the fact that the band with dimers has chains that are twice as long. We then normalise against the sum of the three bands to obtain the relative fraction of chains in each population.

The dynamics of the systems are measured using microrheology. Solutions are made by mixing 8 µl of 1288 linearised plasmid at different concentrations to a final concentration in the range 2ng/µl-500ng/µl with 1 µl of 40 U/ul T4 ligase and 1 µl of T4 ligase reaction buffer. Control solutions are prepared at the same time and in the same manner substituting additional TE for the T4 ligase. The samples are kept at room temperature on a roller for several days. The samples are then spiked with a = 800 nm PVP-coated polystyrene beads, pipetted and sealed onto a slide and imaged using an inverted microscope. We take a 30minute movie and we analyse the movies using a particle tracking algorithm (trackpy [5]) and extract the trajectories and mean squared displacements(MSD) of the tracers ∆r 2 (t) = [r(t + τ )r(t)]

2 . Diffusion coefficients are extracted by fitting to the MSD's via MSD= 2Dt. The viscosity is obtained using the Stokes-Einstein relation [6], η = k B T /(3πDa).

GENERIC CONSIDERATIONS ON THE PROPOSED SMOLUCHOWSKI EQUATION

In this section we explain in detail how we solve the modified Smoluchowski equation reported in the main text. The linear polymers undergo irreversible ligation with rate k 1 (i, j) (i, j=lengths of the filaments reacting) and cyclisation with rate k 0 (l). The concentrations of filaments of length l at time t, ρ(l, t), and of rings, ρ r (l, t), are thus governed by the following equations:

ρ(l, t) = 1 2 ij;i+j=l k 1 (i, j)ρ(i, t)ρ(j, t)+ -ρ(l, t) ∞ i=1 k 1 (l, i)ρ(i, t) -k 0 (l)ρ(l, t) (5) ρr (l, t) = k 0 (l)ρ(l, t) . ( 6
)
Once a chain undergoes cyclisation, it cannot react anymore and thus is removed from the ensemble of linear chains and transferred to the ensemble of rings. The dynamics is also constrained by the requirement that no monomers are created or destroyed and the total mass is conserved:

∞ l=1 l[ρ(l, t) + ρ r (l, t)] = N/V = n ∀t ( 7
)
where N is the total number of monomers and V is the system's volume. Assuming that the reaction takes place on a time scale larger than the Rouse relaxation time, the length-dependence of the annealing rate is [7] k

1 (i, j) = κ 1 (D i + D j )(R i + R j ) (8) = κ 1 l -α i + l -α j l ν i + l ν j , (9)
where κ 1 is a constant coefficient and l i = il 0 is the length of a polymer with degree of polymerisation i, with l 0 the monomer size. The cyclisation rate is taken to be k 0 (l) = κ 0 l µ , where µ = -4ν. Note that we take µ = -4ν instead of the Shimada-Yamakawa exponent µ = -3ν because we are dealing with a non-equilibrium, subdiffusive search process (see below). As we describe in detail later (see "Looping probability and the exponent µ" section), both theory and simulations suggest that the irreversible looping probability scales with the length of the chain as l -4ν and not as l -3ν .

Time scaling of the mean length -concentrated regime

We can solve the Smoluchovski equation in the concentrated limit, where we expect the ring formation to be negligible. We can thus use a mean-field expression for the mean length l and estimate (at scaling level) the rate as

k ∼ DR ∼ l ν-α . ( 10
)
The timescale of chain growth is related to this rate and, in a mean-field picture, we can write dρ(t)/dt = -kρ(t) 2 or k ∼ 1 ρ(t)t (11) according to dimensionality arguments. In these equations, ρ(t) is the number density of chains with mean length l. This also means that

t ∼ 1 ρ(t)k ∼ l -ν+α ρ(t) . ( 12
)
We now point out that the average number of fragments is related to the average length of the fragments as ρ(t) = l 0 ρ(0)/l. This is because mass is conserved so ρ(0)l 0 = i ρ i (t)(il 0 ), or in mean field, ρ(0)l 0 = ρl at any time. Hence we find

t ∼ 1 ρ(t)k ∼ l -ν+α ρ(t) ∼ l -ν+α+1 ρ(0)l 0 ∼ l 1+α-ν . ( 13
)
To find how the mean length grows in time we can just invert this equation to find

l(t) ∼ t 1/(1+α-ν) ∼ t 1/(1-λ) ≡ t γ , (14)
with λ = να. Notice that α = 1 for Rouse dynamics and α = 2 for reptation. Recall that here we are considering the diffusion in 3D of the chain, not its curvilinear diffusion, hence why α = 2 for reptating chains. At the same time, we consider ν = 1/2 for ideal chains and ν = 3/5 for self-avoiding chains. For instance, by plugging in these numbers we obtain l ∼ t γ with γ 0.66 -0.7 for Rouse (so above overlap but not entangled) and γ 0.4 for reptation (so above overlap and entangled).

Time-dependence of the mean length -dilute regime

In the very dilute regime, we can assume that rings and short n-mers should be more favourable. We thus now solve the Smoluchovski equation analytically assuming that only monomers, dimers and monomer rings are present. Let us assume that the concentration of rings, monomers and dimers is denoted by n r , n 1 and n 2 , respectively. The dynamical system that describes the evolution of those populations is the following:

dn r (t) dt = k 0 (1)n 1 (t) (15)
dn 1 (t) dt = -k 1 (1, 1)n 2 1 (t) -k 0 (1)n 1 (t) (16
)

dn 2 (t) dt = 1 2 k 1 (1, 1)n 2 1 (t) (17)
We solve Eq. ( 16) neglecting the second term as n 2 1 1 in the infinite dilution limit:

n 1 (t) = n 1 (0)e -k0(1)t (18)
The concentration of monomer rings is thus

dn r (t) dt = k 0 n 1 (0)e -k0(1)t (19) which yields n r (t) = n 1 (0)(1 -e -k0(l)t ) . (20)
Substituting in Eq. ( 17), we get

dn 2 (t) dt = 1 2 k 11 n 2 1 (t) = 1 2 k 11 n 1 (0)e -k0(1)t 2 (21)
which yields

n 2 (t) = 1 4k 0 k 12 n 2 1 (0) 1 -e -2k0(1)t . (22)
Assuming these three are the only contributions to the system, then the mean length is given by the following relation

l(t) = l 0 n 1 (t) + l 0 n o (t) + 2l 0 n 2 (t) n 1 (t) + n o (t) + n 2 (t) = = l 0 n 1 (t) + n o (t) + 2n 2 (t) n 1 (t) + n o (t) + n 2 (t) . (23)
In denser solutions, where the population is more polydisperse, the Smoluchowski equation cannot be solved analytically and we refer to the previous section for a scaling prediction and in the section "Looping probability and the exponent µ" for a perturbative approach in the limit of small cyclisation. In the next section we explain the numerical method that we implement to solve it.

Time-dependence of the mean length -concentrated regime

To numerically solve the Smoluchowski equation Eq. ( 6) we wrote a code in MATLAB using the explicit Euler scheme. The numerical evaluation of Eq. ( 6) is then iterated to find the best free parameters κ 0 and κ 1 that fit the MD or experimental mean length versus time data. We recall the expression of the ligation rate:

k 1 (i, j) = κ 1 [(il 0 ) ν + (jl 0 ) ν ][(il 0 ) -α + (jl 0 ) -α ]
and of the cyclisation rate k 0 (l) = κ 0 l µ . The fit is done using the nonlinear least squares MATLAB embedded function (i.e. lsqcurvefit). The rate of ring formation κ 0 and the rate of linear chains formation κ 1 are extracted from this fit by considering four independent sets of 10 independent simulations replicas so that we obtain an error on the fits.

Solving the Smoluchowski equation to fit the data from MD simulations consists of two main parts:

1. We create an objective function for the lsqcurvefit (called Obj smoluchowski) that takes as input the array of initial coefficient guess K 0 = (κ 1 , κ 0 ) and the time data array xdata. It returns the average length as a function of time, array ydata. In the objective function:

(a) An array L = n • N 0 is initialised where n = {1, 2, . . . , N m = 200} and N 0 = 174. This represents the set of lengths that can be found in the system (recall that we initialise our MD simulations with 200 chains of 174 beads each). Also, the arrays with the number density of linear and ring chains are initialised as follows, n l0 = (N m /vol, 0, . . . , 0) 1×Nm and n r0 = (0, . . . , 0) 1×Nm since initially all the molecules are linear chains.

(b) for t = {1 to simulation final step time} do call

(n Lnew , n Rnew ) = exEuler smoluchowski (n L , n R , K) func- tion (see point 2 below) (c) update arrays n L = n Lnew and n R = n Rnew .
Calculate the total average length l total as

l total (t) = n Lnew • L + n Rnew • L Nm i n i Lnew + Nm i n i Rnew (d)
exit for loop and parse l total (t) to ydata 2. The exEuler smoluchowski function takes as input the initial number densities of linear and ring chains and the reaction rates n L , n R , K. Based on the given rates K, it outputs final number density arrays n Lnew , n Rnew , after the reactions have taken place. When this function is called, the number density of linear and ring chains of each population are updated according to Eq. ( 6). The monomer, dimer, and so on populations are increased according to the first two terms of Eq. ( 6) while the number of them that is converted into rings is subtracted by the n Lnew and added to the n Rnew array.

In the first two terms of Eq. ( 6) the rate k 1 (i, j) is not a scalar quantity by rather a matrix that follows the relation Eq. ( 9). The extracted coefficient against which the fitting is optimised is the scalar κ 1 . Similarly, for the sink term of Eq. ( 6), the equation k 0 (l) = κ 0 l -4ν is used and the fitting coefficient exported is the scalar κ 0 .

The coefficients K = (κ 1 , κ 0 ) are updated iteratively by the lsqcurvefit algorithm to best fit the data. Once the optimum values are obtained the algorithm terminates.

Looping probability and the exponent µ

In equilibrium the looping probability of a chain of length L is given by the Shimada-Yamakawa formula [8,9]. For L l p (with l p the persistence length) the looping probability decays as P (L) ∼ L µ with µ = -3ν. This looping probability also holds for an irreversible, nonequilibrium scenario if the process is reaction-limited. This is because the chain ends would have the time to explore many conformations and to diffuse the whole volume of the chain, V ∼ L 3ν , before reacting (as it would happen in equilibrium). In a diffusion-limited, irreversible ligation process, one should instead compute the time it takes for an end to diffuse over a certain length distance ξ. The dynamics of the end is described by the Rouse model [10], so that ξ = b(D 1 t/b 2 ) 1/4 , where b is the size of a monomer and D 1 its free diffusion coefficient. Then, setting ξ = R (the size of the chain) one obtains (R/b) 4 = (D 1 t/b 2 ), which implies k 0 ∼ t -1 ∼ R -4 ∼ L -4ν . So considering µ = -4ν effectively takes into account the fact that the chain ends are performing a sub-diffusive search process within the polymer coil, as expected for Rouse dynamics.

We have used MD simulations to check the dependence on the polymer length of our loop formation. To do this we ran short dilute simulations and starting from polymers with different lengths, monitored the number of rings N rings over time and then taken the derivative dN rings /dt in the limit t → 0. As shown in Fig. 1, the rate of ring formation dN rings /dt = k 0 (L) = κ 0 L µ follows a power law which we can fit with an exponent that is close to µ = -4ν, with ν = 0.588 the exponent expected for self-avoiding chains. In light of this evidence, we decided to use µ = -4ν in the main paper.

where we defined x = l/t 1/ (1-λ) . Thus, the number density of polymers that are converted into rings over infinite time is

ρ ∞ r (l) = 2κ 0 (1 -λ)l µ-1-λ ∞ 0 x λ G(x)dx . (32)
Since λ < 1 and assuming the G(x) = O(1) when x → 0, the integral converges in 0. For convergence of this integral at ∞ we also require that the scaling function decays faster than x λ-1 . Assuming this functional form for the distribution of ring lengths at infinite time, we now compute the total average mass transformed into rings at infinite time as

M ∞ rings = ∞ 1 lρ ∞ r (l, t)dl = 2κ 0 (1 -λ) ∞ 0 x -λ G(x)dx ∞ 1 l µ-λ dl . ( 33
)
The convergence of this integral requires that λµ > 1 and in this case we get

M ∞ rings = 2κ 0 1 -λ λ -µ -1 ∞ 0 x -λ G(x)dx . ( 34
)
From this equation we see that the fraction of mass in rings at infinite time M ∞ rings /M 0 converges to a finite value proportional to κ 0 (and hence < 1 at small κ 0 ).

With this calculation we have thus shown that at small enough but non-zero κ 0 , the fraction of mass turning into rings is finite if λµ = να + 4ν > 1 or ν > α/5 which is valid for any type of polymer in the non-entangled (α = 1) regime, and for ideal and self-avoiding polymers in the entangled regime (α = 2). This implies that in these regimes we expect the cyclisation probability to decay fast enough and cannot prevent the runaway of the M 0 -M ∞ rings mass into linear chains that keep growing in time.

Consistently with this, in both MD and MC simulations we never observe the formation of rings larger than 10 initial monomers.

DESCRIPTION OF THE DSMC ALGORITHM

The modified Smoluchowski equation proposed here (6) can only be solved analytically for certain forms of the ligation rate k 1 (i, j) and of the cyclisation rate k 0 (l). Moreover, the computational cost of directly solving the system of ordinary differential equations as described in Sec. increases steeply with the system size. Thus, in order to find solutions to the Smoluchowski equation for larger system sizes, we employ here Direct Simulation Monte Carlo (DSMC) [12][13][14]. DSMC is a powerful stochastic method to solve differential equations such as Eq. ( 6), and which samples the correct ligation kinetics in the limit of large system sizes. The starting point for the Monte Carlo algorithm is an array m of length N , each element i of which contains a number m i which represents the mass/length of the chains i:

m = (m 1 , m 2 , . . . , m N ) .
A value of 0 corresponds to absence of a certain chain. Moreover, to satisfy mass conservation we ensure that N i=1 m i = N (total number of monomers) is true for any time during the simulation. We will also consider an analogous array of length N denoted r (initially empty), where we save the masses of the rings.

For an initially monodisperse condition, we set m 0 = (1, 1, . . . , 1). After the array m is initialized, we run the DSMC simulation, which consists in repeating a large number of times a Monte Carlo step (described below). The execution is terminated when the the system has reached a state in which there is a single linear chain and a number of non-reactive rings.

Description of the Monte Carlo step

In this paragraph, we will describe the single Monte Carlo (MC) step, which is repeated a large number of times during the numerical resolution of Eq. ( 6) performed using the DSMC algorithm. With reference to Eq. ( 6), we define n f ≡ V ∞ i=1 n i (total number of chains). Before the start of the simulation, we give an estimate of the maximum annealing rate k max and of the maximum cyclisation rate k max 0 . The exactness of the algorithm does not depend on this initial choice, however, choosing values that are too far from the actual maximum rates can lead to a reduced efficiency [12].

During every MC step, we either attempt to perform an ligation reaction (with probability p) or a cyclisation one (with probability 1p). The value of p is calculated initially and then updated during the course of the simulation in such a way that the average number densities n(l) satisfy (6). At the beginning of each MC step, p is evaluated as

p -1 = 1 + 2N k max o (n f -1)nk max . (35)
We will show below that this choice also guarantees that the simulation samples the correct number of cyclisation and ligation events per unit volume and unit time as required by Eq. ( 6). We define a waiting time variable that is set to zero at the beginning of the simulation. After each reaction, a waiting time increment is generated. These increments are also chosen in order to guarantee the correct number of ligation and fragmentation reactions per unit time/volume, as detailed below. We can now describe the MC step, during which the following actions are performed: 2. We pick a random number 0 ≤ r ≤ 1 from a uniform distribution. If r ≤ p, we attempt a ligation event:

(a) We pick a pair of elements of the array m, denoted α, β at random. Since there are n f (n f -1) ordered pairs of chains, the probability to pick a specific pair is [n f (n f -1)] -1 . Let the length associated with these elements be m α = i and m β = j.

(b) We evaluate the ligation rate k 1 (i, j) for the two chains. If k 1 (i, j) > k max , we set k max 1 = k 1 (i, j) and return to (1). Otherwise, we continue.

(c) We pick another random number r , and perform the ligation if r ≤ k 1 (i, j)/k max 1

. If ligation is unsuccessful, we return to (1). Otherwise, we continue.

(d) We increment the waiting time by ∆t lig i,j = 2AN n f (n f -1)nkij . Here A is a parameter, the only condition on which is that it must be between 0 and 1, as we will discuss in more detail below.

(e) After incrementing the waiting time, we update f by setting m α = 0 and m β = i + j. We now prove that the definitions of p (Eq. ( 35)), the waiting time increments ∆t lig i,j (for ligation) and ∆t cyc i,k-i (for cyclisation) give a number of ligation and cyclisation events per unit time which is consistent with the Smoluchowski equation Eq. (6). Over a single MC step, the mean number of ligation events involving the ordered pair of filaments (α, β) is

#L α,β ≡ p n f (n f -1) k 1 (m α , m β ) k max 1 . (36)
We note that in the algorithm we consider (m α , m β ) as an ordered pair, and thus in [START_REF] Garcia | [END_REF] we consider the reaction (i, j) → l as distinct from (j, i) → l. The mean number of ligation events involving any two chains with lengths i, j can be obtained by multiplying the above quantity by 2(1δ ij /2)V 2 n i n j . The factor 2(1δ ij /2) takes into account the fact that, as mentioned above, for i = j, there are two ways to perform the ligation, whereas for i = j there is only one. The factor V 2 n i n j is the product of the volume fractions of filaments of lengths i and j. We thus have

2V 2 n i n j 1 - δ ij 2 × p n f (n f -1) k 1 (i, j) k max 1 = V k 1 (i, j)n i n j 1 - δ ij 2 ∆t , (37)
where we have equated the mean number of ligation events involving any two chains with lengths i, j to the value required by the Smoluchowski equation. Recalling that n = N/V , we thus find ∆t = 2pN n f (n f -1)nk max 1 .

(

) 38
Eq. ( 38) relates the time interval ∆t to the probability of ligation. We will now obtain a second equality involving p and ∆t, which will allow us to prove that the expression Eq. ( 35) for p guarantees the correct number of ligation and cyclisation events per unit time.

The mean number of cyclisation events involving chains γ is

#C γ ≡ (1 -p)k 0 (m γ ) n f k max 0 . ( 39
)
To obtain the mean number of cyclisations of a generic l-mer we need to multiply this quantity by V n l , i.e., the volume fraction of filaments of length l. Equating this quantity to the expected number of rings formed in a time interval ∆t we obtain

V n l × (1 -p)k 0 (l) n f k max 0 = k 0 (l)n l V ∆t , (40)
and hence

∆t ≡ 1 -p n f k max o . (41)
By equating the two expressions for ∆t, Eq. ( 38) and Eq. ( 41), we find Eq. [START_REF] Doi | The theory of polymer dynamics[END_REF]. We have thus proven that the latter is the correct expression of p, which gives the correct number of cyclisation and ligation events per unit time and unit volume, as required by the Smoluchowski equation.

Finally, we will prove below that the constants A and 1 -A introduced when calculating the waiting time increments are consistent with Eq. (38) and Eq. (41). To show this, it is sufficient to observe that the total time increment during an MC step is:

∆t = 0≤α<β≤n f -1 #L α,β ∆t lig mα,m β + mγ -1 i=1 #C γ ∆t cyc mγ = 0≤α<β≤n f -1 pk mα,m β n f (n f -1)k max 2AN n f (n f -1)nk mα,m β + mγ -1 i=1 (1 -p)k o (l) nk max o 1 -A n f k o (l) = 2ApN n f (n f -1)nk max + (1 -A)(1 -p) n f k max o ( 42
)
One can see that this equality is consistent with Eq.( 38) and Eq.( 41). We note that the algorithm samples on average the correct kinetics independently of the value of A, as long as 0 ≤ A ≤ 1. Here we take A = 1, meaning that the waiting time increment is calculated only after a successful ligation reaction, but not after a successful cyclisation reaction.

FIG. 1 .

 1 FIG. 1. a. Molecular Dynamics simulations of polymer condensation. b. Snapshot of the simulation box with rings colored in blue and linear chains in grey. c. Pictorial representation of the simulation workflow, tracking the bond list to reconstruct length and topology of condensed polymers. d. Number fraction of linear and ring polymers during the condensation process for different choices of concentration c (averaged across 40 independent replicas). Note that the fraction of rings converge to a finite number at large time. e. Number averaged polymer length for the simulations in d as a function of time, and fitted with Eq. (1). f. Time-resolved gel electrophoresis during ligation of a 2 ng/ul (c/c * = 0.01) solution of 6.5 kb linearised plasmid. The lanes marked with 'e' are treated with exonuclease to remove linear DNA; 'RM' indicates ring monomers, 'LD' linear dimers and 'LM' linear monomers. The last column is λ-HindIII marker to give a reference of DNA fragment lengths. The numbers at the top represent the minutes of ligation treatment. g. Number fraction of polymers in ring and linear (monomer and dimer) topologies. h. Number averaged length calculated from 3 independent gels at 2 ng/µl and associated fit using Eq. (1).

FIG. 2 .

 2 FIG. 2. Runaway transition. a. Dimensionless cyclisation parameter κ = 2κ0/(n0κ1) as a function of c/c * . The scaling κ ∼ (c/c * ) -1 is consistent with κ0 and κ1 being independent on concentration and with this assumption breaking down near c/c * 1. b-d. Snapshots of MD simulations with Nc = 200 monomers: the blue polymers represent rings of any length, grey polymers linear chains.In the network representations each circle i represents a monomer and its colour represents the number of chains connected to i: 0 (grey), 1 (cyan), or 2 (blue). A linear monomer is grey, a ring monomer blue, a linear dimer has two cyan nodes, and a dimer ring would have two blue nodes. e Fraction of monomers in the giant connected component as a function of κ. f. Average polymer length (including both linear and ring) at a large time as a function of κ. g. Average length of only linear chains at a large time as a function of κ. h. Average total number of monomers belonging to rings Mrings = Nrings × lrings(t 1) at large simulation time. Note that some some simulations may not display any rings, in turn bringing the average Mrings below 1. MC = "Monte Carlo"; MD="Molecular Dynamics".

FIG. 3 .

 3 FIG. 3.Rheological consequences of the runaway transition. a. Mean squared displacements (MSDs) of 800 nm PVP-coated polystyrene beads diffusing in solutions of DNA. We compare the MSDs in solutions treated for 1 week either with 40U T4 ligase or buffer (control). b. Viscosity of the ligated solutions η, normalised by the viscosity of the control η0. The inset shows the same plot in log-linear scale to highlight the exponential increase. κ = 1 the GCC displays a change in scaling, growing as GCC ∼ κ -1 ∼ c/c * as κ → 0 (Fig.2e) suggesting that a change in behaviour takes place around κ 1.Since MD simulations are computationally expensive, we employ Direct Simulation Monte Carlo (DSMC)[36- 38] to solve the Smoluchowski equation and reach larger system sizes up to 10 5 chains (see SI). In Fig.2f, we plot the number averaged chain length at an arbitrarily large time when the DSMC code has evolved the system as long as possible and has generated only a single linear chain accompanied by a number of rings. This figure suggests that the linear-dominated regime (κ < 1) displays an average polymer length at large reaction time that scales as l(t 1) ∼ κ -1 ∼ c/c * . At the same time, the number average length of the linear population appears to reach a plateau that scales with the system size (Fig.2g), indicating a runaway transition at κ 1. Additionally the fraction of mass "lost" in forming rings decays as f ring ∼ κ and is thus negligible for small enough κ (Fig.2h).Although our simulations support the notion that a low value of κ will result in linear chains of increasing lengths and vanishing cyclisation rate, they are fundamentally limited to finite-size systems where the cyclisation rate of the largest chain never rigorously goes to 0. To firmly establish the existence of our runaway transition, in the SI we use an asymptotic theory to analytically demonstrate that the mass fraction of linear polymers goes to a finite limit at t → ∞ in a thermodynamic system. We find that the key condition to ensure the runaway transition to exist is that the cyclisation rate k 0 = κ 0 l µ decays strongly enough. More specifically, we require the exponent µ to be µ = -4ν < -4α/5 or ν > α/5. This condition is always met in the Rouse unentangled (α = 1) regime, and also in the entangled (α = 2) regime, provided that the polymers are not fully collapsed (ν = 1/3).Dynamics -To test the consequences of the runaway transition on the dynamics, we perform microrheology experiments (see SI for details) where we track the dynamics of 800 nm PVP-coated polystyrene beads added

FIG. 1 .

 1 FIG. 1. Computing the rate of loop formation as a function of polymer length. (Left) Nrings as a function of time for a system with N0 = 400 chains with M = 87 monomers each. The black line represents the average of the numerical derivatives of the curves in the limit t → 0. (Right) Numerical derivative dNrings/dt showing a power law decay as a function of initial polymer length. The exponent of the power law is found to be -2.6, close to the values -4ν which matches our calculation assuming the irreversible loop formation to be a diffusion-limited and sub-diffusive search process (see text).

3 ., set k max 0 =.

 30 If r > p, we attempt a cyclisation event: (a) We pick a chain γ at random with probability n -1 f . Let m γ = l. (b) We evaluate the cyclisation rate k 0 (l). If k 0 (l) > k max 0 k 0 (l) and return to (1). Otherwise, we continue. (c) We extract another random number 0 ≤ r ≤ 1 from a uniform distribution, and perform cyclisation if r ≤ k 0 (l)/k max 0 If cyclisation is unsuccessful, we return to (1). Otherwise, we continue. (d) We increment the waiting time by ∆t cyc l = 1-A n f ko(l) , with A defined above in step (2). (e) We record the value of l in r and set m γ = 0.

ACKNOWLEDGEMENTS

DM acknowledges the support of the Royal Society via a University Research Fellowship. This project has received support from European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 947918 to DM and No 677532 to ML). The authors acknowledge insightful discussions with Daan Noordermeer and Antonio Valdes. Source codes are available at https: //git.ecdf.ed.ac.uk/taplab/dna-ligation.git.

1. We evaluate the probability of annealing p according to Eq. (35). The explicit form of p, Eq. ( 35), will be discussed in detail below.

Supplementary Information

Maria Panoukidou, 1, * Simon Weir, 1, * Valerio Sorichetti, 2 Yair Gutierrez Fosado, 1 Martin Lenz, 2, 3 and Davide Michieletto 1, 4, †

CALCULATION OF MASS CONVERTED INTO RINGS AT INFINITE TIME

To estimate the amount of mass that is converted into rings at long times, we do a perturbative calculation valid in the limit of small κ. We start from the continuum Smoluchovski equation:

where we define K ≡ k 1 /κ 1 and recall that κ ≡ 2κ 0 /(n 0 κ 1 ) which is thus a scaling function such that K(ai, al) ∼ a λ k 1 (i, l) where λ = να [11].

We now treat κ 0 perturbatively, starting with κ 0 = 0. In this case, there is no mass lost into rings and we can thus write a conservation law

Even for κ non zero, we assume the loss of mass to cyclisation remains finite and of order κ. We will check the self-consistency of this assumption below. Using the mass conservation and Eq. ( 24) we can write the following scaling relations: l 2 ρ = 1, ρt -1 = l 1+λ ρ 2 . We therefore obtain:

which is the same as Eq. ( 14). Note that we must have λ < 1 for the average length of polymers to increase over time. We can also write the density distribution as

which in the limit of long times or large lengths may be written as

where G is a scaling function that only depends on the ratio l/t 1/(1-λ) . We now introduce the ring length distribution ρ r (l, t) and its evolution equation as

Since at time t = 0 there are no rings, we can then write

We can plug in the result we obtained for the distribution of length of linear chains Eq. ( 28) to yield ρ r (l, t → ∞) = 2κ 0 l µ ∞ 0 t -2/(1-λ) G l t 1/(1-λ) dt = 2κ 0 l µ (1λ)l -(1+λ) ∞ 0

x -λ G(x)dx [START_REF] De Gennes | [END_REF]