\section*{Supplemental Material for \\ ``Crosslinker mobility governs fracture behavior of catch-bonded networks''}
\setcounter{equation}{0}
\setcounter{figure}{0}
\setcounter{table}{0}
\setcounter{section}{0}

\renewcommand\thefigure{S\arabic{figure}} 
\noindent

\section{Simulation Methods}
\subsection{Network model}
We perform simulations of 2D triangular networks of $L\times L$ nodes interconnected by one-dimensional Hookean springs with stretching modulus $\mu$ and rest length $l_{0}$. Thus, the Hamiltonian $\mathcal{H}$ of the system is

\begin{equation}
\mathcal{H}=\frac{\mu}{2l_{0}}\sum_{\left\langle ij\right\rangle}\left(l_{ij}-l_{0}\right)^{2}
\end{equation}

\noindent where the sum is taken over all bonds $\left\langle ij\right\rangle$ with $l_{ij}$ the bond length. All bonds are harmonic springs with unit stiffness and unit rest length (so $\mu=1$ and $l_0=1$). We use periodic boundary conditions in the $x$-direction, while nodes on the top and bottom boundaries are fixed. After each deformation step or binding/unbinding event, the energy of the system is minimized by adjusting the positions of the nodes, using the FIRE algorithm~\cite{bitzek2006structural}. The maximum tolerance $F_{RMS}$ corresponds to the system root-mean-square force and it is set to be $10^{-5}$. Since we only consider potential energy, our system is under athermal conditions. Node pairs that have no bond connected do not interact at all. Uniaxial deformation of the network in the $y$-direction leads to a macroscopic  stress in the network $\sigma$. The stress is computed by the virial stress tensor defined as

\begin{equation}
    \sigma_{\alpha\beta}=\frac{1}{A}\sum_{\left\langle ij \right\rangle}f_{ij,\alpha}r_{ij,\beta}\,,
\end{equation}

\noindent where $\alpha$ and $\beta$ are the Cartesian axes. The  sum runs over all the bonded pairs of nodes $\left\langle ij\right\rangle$, $f_{ij,\alpha}$ is the force acting on the node $i$ due to $j$ in the $\alpha-$direction, $r_{ij,\beta}$ is the distance between the two nodes in the $\beta-$direction, and $A$ is the instantaneous area. We evaluate the $\sigma_{yy}$ tensor component. In our simulations,  we maintain a constant uniaxial stress $\sigma=\sigma_{yy}$ by adjusting the strain.

All quantities are expressed in reduced units, i.e. length scales are expressed in terms of $l_0$, forces in terms of $\mu$, and stresses in terms of $\mu/l_0^2$. Time is expressed in units of the rebinding time $1/k_0^b$ (which is the same for slip and catch bonds, see below). 

\subsection{Unbinding and binding bonds}

Structural changes in the network due to binding and unbinding events are simulated by using a kinetic Monte Carlo scheme due to Gillespie~\cite{gillespie1976general}. At every step, we make a list of all possible unbinding and binding events. Based on the list we calculate the total reaction propensity $K_{tot}=\sum_{i}k_{i}$, where $k_{i}$ is the transition rate of a spring corresponding to an unbinding or binding event, as defined below. Which event is the next reaction to occur is drawn randomly, witch each event weighed with its actual reaction rate. The time interval to the next event is drawn also randomly, assuming exponentially distributed events, i.e. $\Delta\tilde{t}=-\textrm{ln}\left(\textrm{rand}\left(0,1\right) \right)/K_{tot}$, with the dimensionless time then given by $t=\tilde{t} k_{0}^{b}$.

All bonds in our networks are transient, reversible bonds with kinetics of binding and unbinding  that depend on the force acting on the bond, $f = \mu\left(l_{ij}-l_{0}\right)/l_{0}$. For the slip bonds, we assume an unbinding rate that decreases exponentially with the force, as predicted by the Bell model~\cite{bell1978models} :

\begin{equation}
k_{u}^{s}\left(f\right)=k_{0}^{s}\exp\left[\frac{f}{f_{e}^{s}}\right]\,,
\end{equation}

\noindent where $k_{0}^{s}$ is the unforced unbinding rate, and $f_{e}^{s}$ is the force where the off rate has fallen with a factor $1/e$. This is related to the activation length $\delta$ of the bond, as $f_{e}^{s}=k_BT/\delta^s$. 
By contrast,  catch bond behavior is described by the two path-way model as~\cite{pereverzev2005two}

\begin{equation}
k_{u}^{c}\left(f\right)=k_{0,1}^{s}\exp\left[\frac{f}{f_{e,1}^{s}}\right]+k_{0,2}^{c}\exp\left[-\frac{f}{f_{e,2}^{c}}\right]\,.
\end{equation}

\noindent The first term, characterized by $k_{0,1}^{s}$ and $f_{e,1}^{s}$, captures the dissociation along a slip-like path at large pulling force $f$, while the second term, characterized by $k_{0,2}^{c}$ and $f_{e,2}^{c}$ describes dissociation along a catch-path decreasing the unbinding rate (and thus, increasing the lifetime) with force for small forces $f$. To mimic catch-bond behavior, it is necessary that $k_{0,1}^{c} > k_{0,2}^{c}$ and $f_{e,1}^{c} < f_{e,2}^{c}$~\cite{pereverzev2005two,van2020chemical}. In Fig. 1(a) in the main text, we show the bond lifetime $t_{u}=1/k_{u}$ as a function of $f$ for both slip and catch bonds. In particular, we decide to fix $\left[k_{0}^{s}, f_{e}^{s}\right] = \left[0.025, 0.6\right]$ for slip bonds, and $\left[k_{0,1}^{s}, f_{e,1}^{s}\right] = \left[0.005,0.6\right]$ and $\left[k_{0,2}^{c}, f_{e,2}^{c}\right] = \left[0.25,0.3\right]$ for catch bonds. Finally, we consider that the binding rate is also load-dependent and the same for both types of bonds, being expressed as:

\begin{equation}
k_{b}\left(f\right)=k_{0}^{b}exp\left[-\frac{f}{f_{e}^{b}}\right]\,.
\end{equation}

\noindent Thus, we assume that the binding rate between a pair $ij$ becomes smaller with increasing distance between the two nodes. Fig. 1(a) in the main text shows the binding lifetime $t_{b}=1/k_{b}$ for the set of values $\left[k_{0}^{b}, f_{e}^{b}\right] = \left[1,0.6\right]$.

\subsection{Immobile and mobile bonds}

Immobile bonds rebind between the same pair of nodes from which they are unbound, providing memory to the system of its initial structure. By contrast, mobile bonds can form in random new locations, involving thus any two nodes. To allow bonds to redistribute, we allow a new bond to form between a pair of nodes that is already linked by one bond, so that double bonds can form.  This implies that, while for immobile bonds the maximum local connectivity is $z^{\left(i \right)}_{max}=6$, for mobile bonds  this is $z^{\left(m \right)}_{max}=2z^{\left(i \right)}_{max}$. In all our simulations, we impose a maximum number of bonds that can be formed $N$, which is  equal to the number of bonds that is present in a network where every pair of neighbouring bonds is connected with exactly one bond. This means that new bonds can only form if there are unbound linkers available.

\newpage

\section{Additional results}

\subsection{Effective viscosity}

Fig.~\ref{fig:Viscosity} shows the effective elongational viscosity computed as $\eta=\sigma/\dot{\epsilon}$ for networks with slip and catch bonds. For networks with slip bonds, Fig.~\ref{fig:Viscosity}(a), we observe that $\eta$
decreases with increasing stress, indicating strain thinning. For immobile bonds $\eta$ is larger than for mobile bonds, which is in line with our finding that the network lifetime is larger for immobile bonds. In particular, for immobile bonds we observe that the viscosity diverges for low stress, which corresponds to the metastable regime.


By contrast, networks with catch bonds manifest strain thickening, Fig.~\ref{fig:Viscosity}(b). Furthermore, unlike what happens for slips bonds, $\eta$ is considerably larger when catch bonds are mobile. This opposite behavior between bonds as a function of mobility strongly suggests that mobility is needed to take advantage of the nature of catch bonds. 

\begin{figure*}[!h]
\includegraphics[width=0.7\linewidth]{FigS1.pdf}
\caption{Effective viscosity $\eta=\sigma/\dot{\epsilon}$, for networks with (a) slip bonds and (b) catch bonds as a function of the stress $\sigma$ and the mobility. \textit{Inset}: Zoom on the viscosity for mobile slip bonds.}
\label{fig:Viscosity}
\end{figure*}

\subsection{Catch bond behaviour on the network}

\begin{figure*}[!h]
\includegraphics[width=0.45\linewidth]{FigS2.pdf}
\caption{Lifetime of a network with irreversible catch bonds  as a function of $\sigma$. Here, the rebinding rate $k_{b}$ is assumed to be zero.}
\label{fig:FigS1}
\end{figure*}

The lifetime of a network with catch bonds without rebinding transitions, i.e. $k_{b}=0$, is shown as a function of stress in Fig.~\ref{fig:FigS1}. We observe how a peak in $\tau_f$ develops at $\sigma\approx0.5$. The presence of this peak resembles the characteristic peak of an individual catch bond. However, it is less pronounced, probably because the stress distribution in the network is heterogeneous. Since  rebinding events promote bond formation at low stress, the echo from the catch bond behavior is shielded in transient networks with reversible bonds.

\subsection{Evolution of the total number of mobile bonds }

As mentioned in the main text, the total number of bonds $n_{b,1}+n_{b,2}$ for both mobile slip and catch bonds remains more or less constant, except close to the moment of fracture, when the network is destabilized and no more bonds can be bound close the cracks to postpone their propagation.

\begin{figure}[!h]
\includegraphics[width=0.9\linewidth]{FigS3.pdf}
\caption{Evolution of the total number of bonds, $n_{b,1}+n_{b,2}$ for (a) mobile slip bonds and (b) mobile catch bonds.}
\label{fig:nb1nb2}
\end{figure}

\subsection{Fraction of final bonds}
In Fig.~\ref{fig:Fractionfinal}(a) we report the fraction of bonds at the moment of fracture, $n_{f}=N_{f}/N$ as a function of stress for immobile bond networks. In particular, for slip bonds, we see a minimum at $\sigma\approx0.5$, i.e. close to the metastability. The presence of this minimum hints at a change in fracture scenario, from  fracture at subcritical stress at low $\sigma$, to fracture dominated by gradual degradation of the network at high $\sigma$. This minimum is shifted to $\sigma\approx0.22$ for immobile catch bonds, leading to a  metatable state at lower stress, and $n_{f}$ is a monotonic function of $\sigma$, indicating that the fracture becomes more brittle as the stress increases. We note furthermore that the fraction of bonds remaining in the network is much smaller for catch than for slip bonds, indicating much more diffuse damage in the case of catch bonds. 

Fig.~\ref{fig:Fractionfinal}(b) shows the fraction of single $n_{f,1}$ and double $n_{f,2}$ bonds for networks with mobile slip and catch bonds. We observe that for intermediate stresses, $n_{f,2}$ is higher for catch bonds, indicating that the redistribution of bonds by the accumulation of double bonds proceeds more efficiently for catch bonds. This is caused by the fact that catch bond networks are stabilized by relocalizing bonds, while slip bond networks are destabilized,

\begin{figure}[!h]
\includegraphics[width=\linewidth]{FigS4.pdf}
\caption{(a) Fraction of the number of bonds at the moment of fracture, $n_{f}=N_{b}\left(t_{f}\right)/N$, for immobile slip and immobile catch bonds as a function of $\sigma$. The vertical dashed line indicates the onset of the metastability. (b) Fraction of final single $n_{f,1}$ and double $n_{f,2}$ bonds for networks with mobile slip and mobile catch bonds.}
\label{fig:Fractionfinal}
\end{figure}

\subsection{Unbinding and binding probability distribution}
We record the bond length associated to each binding/unbinding event during the simulation, and build the probability distribution that a bond will unbind $P_{u}$ or bind $P_{b}$ with length $l$. This information is shown in Fig.~\ref{fig:Prob} for the four different types of networks we have simulated here, for different $\sigma$. While the mobility of bonds strongly influences the probability distributions for slip bonds, for catch bonds they are quite similar. Indeed, we observe that both $P_{b}$ and $P_{u}$ for immobile slip bonds develop a double peak; the first placed at the bond rest length, whereas the second one shifts to higher bond length with $\sigma$. Because immobility limits the places where bonds can be formed, it is expected that slip bonds are forced to bind with $l_{ij}\gg l_{0}$. However, these strongly stretched bonds will also break rapidly due to their slip behaviour. When slip bonds are mobile, the second peak in $P_{b}$ disappears, and the distribution develops a long tail. Now, slip bonds can form wherever they want, and, since their lifetime is maximum at $l_{ij}\rightarrow l_{0}$ they will preferentially want to accumulate  in regions with lower local stress, as we have also observed in Fig. 3(e) in the main text where the distribution of double bonds around cracks is discussed. This has a direct impact on $P_{u}$, where the second peak at large bond length is reduced. As we have indicated above, the binding probability distribution for immobile and mobile catch bonds are similar.
The unbinding probabilities at higher force are much smaller for catch bonds than for slip bonds, clearly showing that catch bonds are stabilized in high stress regions. 

\begin{figure*}[!h]
\includegraphics[width=0.9\linewidth]{FigS5.pdf}
\caption{Unbinding $P_{u}\left(l\right)$ and binding $P_{b}\left(l\right)$ probability distributions for different stress values $\sigma$. (a) Immobile slip bonds, (b) immobile catch bonds, (c) mobile slip bonds, and (d) mobile catch bonds.}
\label{fig:Prob}
\end{figure*}

\newpage

\subsection{Crack identification}

To identify crack and study crack propagation, we define on the undeformed network, i.e. the fully connected triangular lattice, centroid points of each triangular cell as shown in Fig.~\ref{fig:CrackId}(a). Unbinding and binding events are mapped on the undeformed network, schematically shown in Fig.~\ref{fig:CrackId}(b). Then, two centroid points are linked defining a crack when, during the deformation, the bond that splits them is broken (see Fig.~\ref{fig:CrackId}(c)). Furthermore, we define the crack surface as the lattice nodes or bonds around each crack, as we represent in Fig.~\ref{fig:CrackId}(d). Thus, we compute the crack size distribution $n_{c}$, which we discuss in the main text, as well as the crack surface size distribution $n_{s}$ shown  in Fig.~\ref{fig:SurfaceCrack}(a), exhibiting both the same trend. We also compute the distribution of double bonds around the cracks $P(N_{b,2})$ for different $\sigma$, represented in Fig. 3(e) in the main text.

\begin{figure}[!h]
\includegraphics[width=0.8\linewidth]{FigS6.pdf}
\caption{Sketch representing the crack identification procedure. (a) Undeformed triangular network, where centroid points of the triangles are represented. (b) Broken bonds are mapped on the undeformed network. (c) Once broken bonds are localized, centroid points are linked, defining a crack. (d) Surface nodes and surface bonds are highlighted around the crack defined on the undeformed network.}
\label{fig:CrackId}
\end{figure}

\subsection{Network evolution}
\subsubsection{Surface size distribution}
The crack surface area  distribution $P\left(s\right)$, shown in Fig.~\ref{fig:SurfaceCrack}(a) exhibits the same behaviour as the the crack size distribution represented in Fig. 3 in the main text. Indeed, $P\left(s\right)$ shows a power-law dependence at large crack sizes with an exponent $\alpha\approx-2.2$ just before the network fractures for networks with immobile slip, immobile catch and mobile catch bonds. This value is consistent with the random percolation value~\cite{torquato2002random,sciortino2005one}. In addition, both $P\left(c\right)$ and $P\left(s\right)$ evolve in the same way as a function of $t$. 

\begin{figure}[!h]
\includegraphics[width=\linewidth]{FigS7.pdf}
\caption{Evolution of the surface size distribution $P\left(s\right)$ as a function of the network structure at stress $\sigma=0.36$. In each panel, the full line represents the function $P\left(s\right)\propto s^{-2.2}$. Here, $P\left(s\right)$ is represented at different $t$ in log scale. Only the distributions at early stage ($\left< t\right> \sim t_{ini}$) and close to fracture ($\left< t \right> \sim t_{f}$) are indicated by symbols, to highlight the structural changes.}
\label{fig:SurfaceCrack}
\end{figure}


\subsubsection{Network evolution at low stress}

In Fig.~\ref{fig:Cluster010} we show $P_{c}$ at $\sigma=0.10$ for networks with immobile slip bonds, immobile catch bonds and mobile catch bonds (note that networks with immobile slip bonds do not break at this stress). The distributions look similar to those at higher stress in Fig. 3 in the main text, albeit it that the distributions for catch bonds are even wider at this lower stress.

\begin{figure}[!h]
\includegraphics[width=0.8\linewidth]{FigS8.pdf}
\caption{Evolution of the crack size distribution $P\left(c\right)$ as a function of the network structure at stress $\sigma=0.10$, for immobile catch bonds, mobile slip bonds, and mobile catch bonds. In each panel, the full line represents the function $P\left(c\right)\propto c^{-2.2}$. Here, $P\left(c\right)$ is represented at different $t$ in log scale, although only at early stage ($\left< t\right> \sim t_{ini}$) and close to fracture ($\left< t \right> \sim t_{f}$) are indicated by symbols to highlight the structural changes.}
\label{fig:Cluster010}
\end{figure}

\subsubsection{Double bond distribution}

In Fig.~\ref{fig:DoubleBondsStress} we show the distribution of double bonds $P\left(N_{b,2}\right)$ computed on the crack surface for different $\sigma$ for networks with mobile slip and mobile catch bonds. In particular, for mobile slip bonds we can see that the distribution decays faster by increasing $\sigma$, revealing that bonds are less stable due to the high local stress on the crack surface. By contrast, the effect of stress on $P\left(N_{b,2}\right)$ is smaller for networks with catch bonds. We fit the tail of these distributions by considering $P\left(N_{b,2}\right)\propto N_{b,2}^{-\alpha}$. Te inset of Fig.~\ref{fig:DoubleBondsStress} clearly shows that the accumulation of catch bonds is more pronounced than for slip bonds. Thus, crack propagation is hindered more efficiently, postponing the fracture because there is an accumulation of bonds around the cracks. 


\begin{figure}[!h]
\includegraphics[width=0.8\linewidth]{FigS9.pdf}
\caption{Distribution of double bonds $P\left(N_{b,2}\right)$ computed on the crack surface as a function of $\sigma$ close to the moment of fracture $t_{f}$, for networks with mobile slip bonds (blue) and mobile catch bonds (red). The right arrow highlights the stress direction. {\it Inset:} Exponent corresponding to the fitting $P\left(N_{b,2}\right)\propto N_{b,2}^{-\alpha}$ made to the tail of the distribution.}
\label{fig:DoubleBondsStress}
\end{figure}

\newpage

\subsubsection{Damage visualization}

The evolution discussed by computing $P\left(c\right)$ at $\sigma=0.40$ is represented by snapshots taken at different times for networks with immobile bonds in Fig.~\ref{fig:Evol_imm} and for networks with mobile bonds in Fig.~\ref{fig:Evol_mob}, where single bonds are shown in orange, and double slip and catch bonds are shown in blue and purple, respectively. Likewise, we also visualize the damage propagation on the network once the fracture has taken place, i.e. $t=t_{f}$.

\begin{figure*}[!h]
\includegraphics[width=0.8\linewidth]{FigS10.pdf}
\caption{Network evolution as a function of time $t$ with immobile bonds, for $\sigma=0.40$.}
\label{fig:Evol_imm}
\end{figure*}

\begin{figure*}[!h]
\includegraphics[width=0.8\linewidth]{FigS11.pdf}
\caption{Network evolution as a function of time $t$ with mobile bonds, for $\sigma=0.40$. Orange indicates single bonds, whereas blue and purple indicate double slip and catch bonds, respectively.}
\label{fig:Evol_mob}
\end{figure*}

\begin{figure*}[!h]
\includegraphics[width=\linewidth]{FigS12.pdf}
\caption{Visualization of the damage once the fracture has taken place, i.e. $t=t_{f}$, mapped on the undeformed network at $\sigma=0.10$ (top) and $\sigma=0.40$ (bottom). Here, broken bonds belonging to the percolated crack are highlighted in red color, whereas blue color indicates individual broken bonds. Note that the network with immobile slip bonds does not fracture for $\sigma=0.10$.}
\label{fig:Evol_imm2}
\end{figure*}

\clearpage
\bibliographystyle{unsrt}
