# DiffMoog: a Differentiable Modular Synthesizer for Sound Matching

###### Abstract

This paper presents DiffMoog - a differentiable modular synthesizer with a
comprehensive set of modules typically found in commercial instruments. Being
differentiable, it allows integration into neural networks, enabling automated
_sound matching_ , to replicate a given audio input. Notably, DiffMoog
facilitates modulation capabilities (FM/AM), low-frequency oscillators (LFOs),
filters, envelope shapers, and the ability for users to create custom signal
chains. We introduce an open-source platform that comprises DiffMoog and an
end-to-end sound matching framework. This framework utilizes a novel _signal-
chain loss_ and an encoder network that self-programs its outputs to predict
DiffMoogs parameters based on the user-defined modular architecture. Moreover,
we provide insights and lessons learned towards sound matching using
differentiable synthesis. Combining robust sound capabilities with a holistic
platform, DiffMoog stands as a premier asset for expediting research in audio
synthesis and machine learning. Our code is released at:
https://github.com/aisynth/diffmoog.

Index Terms—  differentiable synthesis, sound matching

##  1 Introduction

Synthesizers are electronic musical instruments capable of generating a vast
spectrum of sounds, from simple tonalities to complex auditory textures,
central to many music genres. Composed of modules like oscillators, filters,
low frequency oscillators, ADSR envelopes and modulators, they are controlled
by an interface to fine-tune the sound output [1]. However, the multitude of
interacting parameters makes sound design complex, often demanding deep
expertise in sound synthesis and iterative testing.

Neural networks are increasingly used in sound design to replicate input
sounds. While initial methods optimized a loss over synthesizer parameters [2,
3, 4], it may be beneficial to optimize over the sound directly, since the
ultimate goal in sound matching is a high fidelity reproduction. Furthermore,
replicating unlabeled out-of-domain sounds not generated by the synth at hand
requires unsupervised learning. The implementation of traditional synthesizers
does not support automatic differentiation, which is vital for direct sound
comparison and optimization via backpropagation of gradients. To address this,
several works proposed differentiable synthesizers [5, 6, 7]. However they
either presented complex [5] or overly simplistic [6] models that deviate from
real world synthesizers. For example, both synthesizers from [5, 6] lack
modularity and conventional sound modules (e.g., LFOs, FM/AM modulators),
leaving a gap for practical applications.

We present DiffMoog - a differentiable synthesizer with a modular
architecture. It incorporates modules commonly found in commercial
instruments, following an explainable design familiar to those versed in sound
synthesis. The modularity enables both the creation of custom signal chains
and the ability to isolate modules for research purposes. DiffMoog pairs with
an end-to-end platform, facilitating its integration into an automatic sound
matching system, along a newly crafted ‘signal-chain loss’ aimed at guiding
the optimization process. The platform supports dataset creation, experiment
configuration, model training, and offers comprehensive logging and analysis
tools. By releasing DiffMoog as open-source, we aim to propel AI-guided sound
synthesis forward. We discuss its implementation, evaluate its potential for
sound matching, and share key insights and lessons gained throughout the
experimentation process. We note that in the context of sound matching, the
sub-task of frequency estimation through gradient descent techniques via
minimizing spectrogram-based losses is an intrinsic challenge that remains
open [8], as we discovered through our own experimentation. Previous studies
[6, 7, 5] have relied on alleviating assumptions about the data or employed
pitch estimation algorithms like CREPE [9]. However, we did not employ such
techniques in our optimization process.

Our main contributions include: (1) The open-sourced DiffMoog synthesizer and
sound matching platform, aiming to provide an easy gateway to conduct research
in the field of AI sound synthesis and sound matching. (2) The introduction of
a novel signal-chain loss. (3) Lessons and insights learned from optimizing
DiffMoog; and (4) Demonstrating the superiority of the Wasserstein loss in
frequency estimations.

##  2 Related Work

Non-Differentiable Synthesizer Sound Matching: With the rise of machine
learning methods, works in sound matching, also known as _parameters
inference_[3], utilized supervised datasets of sound samples and their
parameters, derived from typical non-differentiable synthesizers. Using this
data, neural networks were trained to predict sound parameters, a method seen
in commercial VST instruments [4, 10, 3]. Sound matching was also explored on
custom synthesizers [11, 2]. While supervised methods optimize over
synthesizer parameters, the ultimate goal of sound matching is obtaining a
high-fidelity sound reconstruction, suggesting a direct audio optimization.
These methods are also bound to synthesizer-specific data. DiffMoog bypasses
these restrictions, being differentiable, allowing direct input-output sound
comparisons beyond parameter reliance, broadening training and evaluation
scope to include out-of-domain sounds. Noteworthy, direct sound optimization
can also be achieved using genetic algorithms (GA) [4, 12]. However, GAs are
inefficient, necessitating a vast number of synthesizer renders for a single
match.

Neural Network Based Synthesizers: Advancements in deep learning gave rise to
a suite of neural network based synthesizers [13, 14, 15, 16, 17, 18]. Unlike
traditional methods relying on the superposition of sinusoids or the routing
of audio submodules (e.g., Additive, Subtractive, FM, Wavetable synthesis and
their variants), these synthesizers generate sound directly as the output of a
neural network. However, they lack the control and explainability that the
classic synthesizers offer. Other works explored the power of generative
adverserial networks . GANSynth [15, 19] and diffusion models [16, 20].
Setting itself apart, DiffMoog uses traditional synthesizer modules, granting
full user control and explainable behavior.

Differentiable DSP: Differentiable digital signal processing (DDSP) [5]
integrates signal processing modules as differential operations into neural
networks, allowing backpropagation. It uses a flexible additive synthesizer
that can produce complex, realistic sounds. Additive synthesis is based on the
Fourier theorem of stacking sinusoidal waves to construct a complex sound. A
DDSP-inspired wavetable synthesizer was introduced in [21], where arbitrary
waveshapes are stacked to generate the final sound. Both methods have limited
control due to the large number of parameters with unpredictable effect on the
final sound, though DDSP does allow manipulation of fundamental frequency and
sound volume. Differentiable methods had also been employed in audio effects
applications [22, 23]. Recent advancements include a differentiable mixing
console for automatic multitrack mixing [24], transferring audio effects and
production style between recordings through differentiable effects [25], and
automating DJ transitions with differentiable audio effects [26].

Differentiable Synthesizers: Differentiable synthesizers aim to model
commercial synthesizer characteristics, often using subtractive synthesis, in
which an initial harmonically rich oscillator is further manipulated to modify
its frequency content. A basic synthesizer, built upon two additive
oscillators and a filter is used in [6]. Recently, the same authors enhanced
their synthesizer to include ADSR envelopes and effects [27]. Yet, they lack
certain complexities found in standard instruments, like LFOs, FM modules, and
flexible routing options. FM synthesis as advised by Chowning [28], uses sine
wave oscillators that are interconnected such that they serve as carriers and
modulators, controlled by ADSR envelopes that can create complex harmonics and
timbres evolving in time. Caspe et al. [7], manually implemented a few
differential patch configurations from the well-known DX7 FM synthesizer [29].
Our proposed DiffMoog, to the best of our knowledge, is the first and most
comprehensive modular differentiable synthesizer of its kind. It integrates
both FM and subtractive synthesis techniques, offering unparalleled
capabilities for sound generation and customization within the domain of
differentiable synthesis.

##  3 The DIFFMOOG Synthesizer and Sound Matching Platform

For clarity, we have aligned the nomenclature in this paper with the terms in
our codebase. The code is organized as follows:
(1) Core implementation of the synthesizer - src/synth.
(2) Functionality of the sound matching system - src/model.
(3) Operations related to dataset management - src/dataset.
(4) Configurations for training, synthesis, and loss functions - configs
directory.

###  3.1 DiffMoog

DiffMoog is a modular synthesizer, inspired by Robert Moog’s seminal Moog
Synthesizer of 1964, and stands as a differentiable synthesizer, optimized for
gradient-based computations. DiffMoog integrates the qualities of subtractive
modular synthesis with rich-harmonic oscillators, while also offering the
capabilities of FM and additive synthesis. Through its versatile array of
modules—oscillators, LFOs, filters, modulators, and shapers — DiffMoog is
capable of producing a wide range of soundscapes surpassing previous
differentiable synthesizers.

The architecture is built upon a 2D _matrix_ of _cells_. These cells act as
containers for the synthesizer _modules_ and their associated _parameters_ ,
as defined in the src/synth/synth_architectur-
e.py file. This matrix is structured by channels (rows) and layers (columns).
Each cell, uniquely identified by its channel and layer indices, can host a
single module. Signal flow within this structure is managed by: (1)
_connections_ which are assets of the cells: each cell is allowed to have
multiple input and output connections, which dictate how signals traverse the
matrix, and (2) using activation parameters to control module behavior,
ranging from a simple on/off switch to bypassing an internal procedure. For
example, one could opt to use just one oscillator where two are present or
neutralize FM in a modulator. In order to simplify the signal flow and
precisely define the sound generation process described later, the
architecture imposes a rule on cell connections: a cell can only receive input
from cells in preceding layers and send outputs to subsequent layers.

Fig. 1: The DiffMoog synth with an arbitrary _chain_. Shown: _matrix_ ,
_cells_ , _modules_ , _connections_ (arrows) and _parameters_ (‘p’). Black
arrows are _fixed connections_ , gray arrows are _optional connections_. Empty
cell outputs the 0 signal.

DiffMoog allows users to craft signal chains using predefined modules and
connections, collectively termed _chains_ , detailed in
src/synth/synth_chains.py. This feature enables versatile synthesizer
configurations, but also permits module isolation for research needs.
Utilizing a particular chain, one can generate random datasets. Within a
chain, connections can be (1) _fixed_ , always present, or (2) _optional_ ,
which, during dataset creation, are decided upon randomly for each sound
instance in conjunction with the module’s activation parameters to finalize
the signal flow. Noteworthy, sound modules may possess optional inputs, and
these must correspond with connections leading to their respective cells. An
arbitrary chain of the DiffMoog synthesizer is depicted in Figure 1.

After populating the synthesizer matrix with modules and defining the cells
interconnections (i.e. setting up a chain), the signal chain and system state
are solidified. Sound generation proceeds layer-by-layer, from the lowest to
the highest order. In each layer, cells compute sound independently,
processing their (optional) input signals via their module logic and
parameters. The signals from the cells then serve as input for cells in
subsequent layers. Upon completion, the final layer’s outputs across all
channels are averaged to yield the final sound output (See
src/synth/synth_architecture.py::SynthModula-
r.generate_signal()).

The synthesizer modules, are implemented in src/synth/s-
ynth_modules.py. The Oscillator generates sine, square, and sawtooth signals
spanning the human hearing range (20Hz-20KHz), employing closed-form
expressions instead of the additive synthesis used in previous works like [6].
The Low-Frequency Oscillator (LFO) acts as a sub-audible control signal,
enriching temporal variations in sound, especially in the FM and Tremolo
contexts. Furthermore, the system integrates FM Oscillator for frequency
modulation, and a Filter typical for subtractive synthesis. Control of
amplitude and filter characteristics over time is facilitated by the ADSR
modules, while the Mix module blends input sounds, and Tremolo modulates
amplitude for a pulsing effects.

(a) AM Square wave using tremolo effect

(b) FM square modulated by a sawtooth wave

Fig. 2: Spectrograms sounds synthesized with DiffMoog. While typical, FM/AM
sounds cannot be synthesized by prior differentiable synths.

Spectrogram plots of several sound examples generated by the DiffMoog
synthesizer are shown in Figure 2. The corresponding audio examples are
available for listening on the DiffMoog’s GitHub repository.

###  3.2 End-To-End Sound Matching Platform

To train a model, an audio input is fed into a neural network, which predicts
parameters for modules in a specified DiffMoog chain. In turn, the predicted
parameters are used to synthesize a sound output that is optimized to
replicate the input. This process is guided by minimizing both spectral and
parameter losses. Since the parameter loss is optional, the system allows
unsupervised training of unlabeled sounds, including out-of-domain sounds not
made by DiffMoog. At inference time, given a new audio input, the neural
network outputs parameter predictions for the corresponding DiffMoog chain
that can be used to synthesize the predicted sound. A system diagram is shown
in Fig. 3.

The encoder network features a neural backbone (e.g., ResNet, GRU) that
produces a latent vector, complemented by a dynamically configured Multi-Layer
Perceptron (MLP) head for each parameter, tailored to the selected synth
chain. The encoder architecture is illustrated in Fig. 4. Notably, connections
are not predicted as they are assets of the cells. Instead, the signal flow is
governed exclusively by classified activation parameters. During training,
connections deemed as ’optional’ within the chain are treated as ’fixed’.
Continuous parameters are normalized using the sigmoid function, while the
‘Gumbel-Softmax’ [30] trick is employed for sampling from categorical
parameters. The reader is referred to src/model/model.py in our repository for
the exact code implementation

The loss function we employed is a combined loss consisting of two terms: a
parameters loss and a newly proposed spectral loss we name _signal-chain
loss_. The parameters Loss, denoted as $\mathcal{L}_{\mathrm{p}}$, is the sum
of synthesizer parameters differences between predicted and original values.
Regression parameters employ L1 or L2 loss, while categorical parameters use
cross-entropy loss. The overall parameters loss is formulated as follows:

$\small\mathcal{L}_{\mathrm{p}}=\sum_{\mathrm{n}\in\mathcal{N}}L_{\text{reg}}(p_{n},\hat{p}_{n})+\sum_{\mathrm{m}\in\mathcal{M}}L_{\text{cat}}(c_{m},\hat{c}_{m})$
(1)

where regression and categorical parameters belong to the sets $\mathcal{N}$
and $\mathcal{M}$ respectively, $r_{n}$ and $c_{m}$ represent the genuine
values, while $\hat{r}_{n}$ and $\hat{c}_{m}$ stand for the predicted values.
The loss functions for these parameter types are symbolized by
$L{{}_{\text{reg}}}$ and $L{{}_{\text{cat}}}$, respectively. The signal-chain
loss measures the difference between the ground truth and predicted audio
signals. However, unlike prior works, it evaluates the signal at all stages
within the synthesizer signal-chain, not only at the final output. Given that
the synthesizer’s output can be seen as a function composition (e.g.,
$y=f(g(x))$ where $y$ is the output signal, $f$ and $g$ are sound modules and
$x$ is an input signal), it aims to guide the optimization by improving early
stage predictions. The signal-chain loss $\mathcal{L}_{\mathrm{SC}}$ is
defined using Eq. 2 and 3 below :

$\small\mathcal{L}_{\mathrm{SC}}=\sum_{\mathrm{i}\in\mathcal{I}}\sum_{\mathrm{j}\in\mathcal{J}}\sum_{\mathrm{k}\in\mathcal{K}}L_{ijk}$
(2)
$L_{ijk}=\left\lVert(F_{k}(S_{ij}(x))-F_{k}(S_{ij}(\hat{x})))\right\rVert_{p}$
(3)

where indices $\mathrm{i}$, $\mathrm{j}$, and $\mathrm{k}$ indicating cell,
FFT window size, and processing type, respectively. The processing function
$F_{\mathrm{k}}$ can be identity, log, or cumulative sum along time/frequency
axes (Wasserstein). The $p$ norm is either 1 or 2, and $S$ denotes a
spectrogram or mel-spectrogram transform. Importantly, the sets
$\mathcal{I},\mathcal{J},\mathcal{K}$ are configurable, allowing various loss
configurations e.g. multiresolution comparison, using traditional output-only
optimization, etc. In total, the combined loss is defined as:

$\small
L_{\mathrm{total}}=\mathcal{L}_{\mathrm{p}}+\beta\cdot\mathcal{L}_{\mathrm{SC}}$
(4)

with $\beta$ as a weighting factor. Please refer to src/model/loss for code
implementation.

Fig. 3: The end-to-end sound matching system diagram.

Fig. 4: The neural network architecture with dynamically allocated MLP heads

##  4 Sound Matching Evaluation

Evaluating DiffMoog for sound matching has been a challenging journey,
encountering non-convergence in many experiments. Despite this, we
successfully established several findings and techniques, presenting a
comprehensive synthesizer and platform as a progressive step in the research
domain. Our work spanned various synthesizer chains, loss configurations, and
neural architectures. While the complexities led to a focus on key findings
rather than exhaustive experimental details, we believe this approach
illuminates the most promising aspects of our research, fostering further
exploration in this field.

Generally, we followed training procedures, hyperparameters and configurations
similar to those outlined in [6]: (1)‘P-loss’: Training with parameters only;
(2)‘Synth’: Training with parameters loss and gradually shifting to spectral
(i.e signal-chain or some other configuration) loss over in-domain data;
(3)‘Real’: Same as (2), but continuing training with spectral loss over out-
of-domain data (Nsynth [13]).

Using the signal-chain loss solely failed systematically, However, when
applied to a relatively basic synthesizer chain (Fig. 5), training with the
spectral signal-chain loss after using the parameter loss hinted superior
performance over the sole usage of parameters loss on out-of-domain data,
which is on par with previous studies [6, 27]. Figure 6 presents an evaluation
example using a chain which forms a synth similar to the one presented in [6],
with a mix of a square and sawtooth oscillators followed by an ADSR amplitude
shaper and a filter.

We also report that more complex chains utilizing FM modulations refused to
converge, adding up to the already hard task of frequency estimation using
spectral loss [8]. The back-propagated gradient from the spectral distance
(Eq. 3) for frequency components and the FM modulation index is very abrupt.
This phenomenon is conveyed in Fig. 7.

In another experiment, we examined the effectiveness of various spectral
processing functions for frequency estimation, adopting the same experimental
setup as described in [8]. A deliberately perturbed signal, more distant from
a target than a predicted signal, is expected to have a gradient oriented
towards the target. This assumption is compared for different loss
configurations. As Table1 shows, using the Wasserstein distance on the time
axis significantly enhanced the accuracy of frequency estimation for both
square and sawtooth waveforms. This outcome defied our initial hypothesis of
seeing more benefits when applying Wasserstein distance on the frequency axis.
Notably, this approach varies from earlier studies that relied on spectrograms
and their log representations for spectral loss

To conclude, our study indicates that differentiable synthesizers equipped
with spectral loss optimization may indeed facilitate sound matching. Yet,
achieving high precision in imitating typical sounds remains a formidable
challenge. Notably, employing the Wasserstein distance could potentially
mitigate the gradient issues encountered in frequency estimation using
spectral loss. We anticipate that our platform will catalyze further research
in this captivating domain. Moving forward, we propose the exploration of
refined audio loss functions, optimization strategies, and alternative neural
network architectures to surmount this hurdle.

Fig. 5: The chain used for the experiment in Fig. 6, with a sawtooth
oscillator, square oscillator, Amplitude ADSR and a Lowpass Filter.

Fig. 6: Training procedures evaluation on NSynth. ’Synth’ uses parameter loss
until step 2k, then gradually introduces spectral loss (trained on in-domain
data) until step 6k. ’Real’ mirrors ’Synth’ but transitions to out-of-domain
data at step 6k. Log-Spectral Distance:
$LSD=\left\lVert(log(S(x))-log(S(\hat{x}))),\right\rVert_{F}$, $F$ is the
Frobenius norm. The ’Real’ procedure shows the efficacy of using out-of-domain
data for training, enhancing real-world sound reproductions.

(a) Loss vs. Oscillator Frequency

(b) Loss vs. Modulation Index

Fig. 7: Loss surface for frequency and modulation index, illustrating highly non-convex behavior with many local minima, which pose challenges for optimization. Other synth parameters are fixed to their ground truth. Loss conf. \ Perturbed dist. | f ± $\varepsilon$ | f ± 300 cents | f ± 600 cents
---|---|---|---
Square waves
Spectrogram, $I$ | 0.5 | 0.676 | 0.7
Mel, $I$ | 0.501 | 0.498 | 0.475
Spectrogram, Wasserstein time | 0.48 | 0.733 | 0.748
Spectrogram, Wasserstein freq. | 0.528 | 0.605 | 0.623
Mel, Wasserstein time | 0.488 | 0.466 | 0.452
Mel, Wasserstein frequency | 0.494 | 0.696 | 0.645
Sawtooth wave
Spectrogram, $I$ | 0.47 | 0.642 | 0.701
Mel, $I$ | 0.52 | 0.506 | 0.463
Spectrogram, Wasserstein time | 0.48 | 0.712 | 0.715
Spectrogram, Wasserstein freq. | 0.53 | 0.571 | 0.543
Mel, Wasserstein time | 0.469 | 0.487 | 0.467
Mel, Wasserstein frequency | 0.505 | 0.624 | 0.619
Table 1: Comparison of different Loss configurations (modified $S$,
$F_{\mathrm{l}}$, in Eq. 3) and waveshapes across various distances of the
pertrubed signal. In a well behaved distance, the prediction should be closer
to the target than the perturbation is. This condition gives us a 0/1 error.
Repeated and averaged over 1000 trials, high accuracies indicate that
gradients usually point in the right direction. f is the predicted wave
frequency, $\varepsilon$ represents the local gradient, $I$ is the identity
function.

## References

  * [1] M. Russ,  Sound synthesis and sampling,  Taylor & Francis, 2004.
  * [2] O. Barkan, D.Tsiris, N. Koenigstein, and O. Kats,  “Inversynth: Deep estimation of synthesizer parameter configurations from audio signals,”  IEEE/ACM Trans. on Audio, Speech, and Language Processing, vol. 27, no. 11, pp. 2385–2396, 2019.
  * [3] P. Esling, N. Masuda, A. Bardet, R. Despres, and A. Chemla-Romeu-Santos,  “Universal audio synthesizer control with normalizing flows,”  in International Conference on Digital Audio Effects (DaFX 2019), Birmingham, United Kingdom, Sept. 2019,  DaFX 2019.
  * [4] M. J. Yee-King, L. Fedden, and M. D’Inverno,  “Automatic programming of vst sound synthesizers using deep networks and other techniques,”  IEEE Trans. Emerg. Topics Comput. Intell., vol. 2, no. 2, pp. 150–159, April 2018.
  * [5] Jesse Engel, Lamtharn (Hanoi) Hantrakul, Chenjie Gu, and Adam Roberts,  “Ddsp: Differentiable digital signal processing,”  in International Conference on Learning Representations, 2020.
  * [6] Naotake Masuda and Daisuke Saito,  “Synthesizer sound matching with differentiable dsp.,”  in ISMIR, 2021, pp. 428–434.
  * [7] Franco Caspe, Andrew McPherson, and Mark Sandler,  “Ddx7: Differentiable fm synthesis of musical instrument sounds,”  arXiv preprint arXiv:2208.06169, 2022.
  * [8] J. Turian and M. Henry,  “I’m sorry for your loss: spectrally-based audio distances are bad at pitch,”  arXiv preprint, vol. arXiv:2012.04572, December 2020,  Published in I Can’t Believe It’s Not Better! (ICBINB) NeurIPS 2020 Workshop.
  * [9] Jong Wook Kim, Justin Salamon, Peter Li, and Juan Pablo Bello,  “CREPE: A convolutional representation for pitch estimation,”  in 2018 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP 2018, Calgary, AB, Canada, April 2018, IEEE, pp. 161–165.
  * [10] K. Itoyama and H. G. Okuno,  “Parameter estimation of virtual musical instrument synthesizers,”  in Proc. 40th Int. Comput. Music Conf., 2014, pp. 1426–1431.
  * [11] Oren Barkan and David Tsiris,  “Deep synthesizer parameter estimation,”  in ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2019, pp. 3887–3891.
  * [12] Kıvanç Tatar, Matthieu Macret, and Philippe Pasquier,  “Automatic synthesizer preset generation with presetgen,”  Journal of New Music Research, vol. 45, no. 2, pp. 124–144, 2016\.
  * [13] J. Engel, C. Resnick, A. Roberts, S. Dieleman, D. Eck, K. Simonyan, and M. Norouzi,  “Neural audio synthesis of musical notes with wavenet autoencoders,”  in Proc. of the 34th Int. Conf. on Machine Learning, 2017, pp. 1068–1077.
  * [14] A. Défossez, N. Zeghidour, N. Usunier, L. Bottou, and F. Bach,  “SING: Symbol-to-instrument neural generator,”  in Advances in Neural Information Processing Systems, 2018, pp. 9041–9051.
  * [15] Jesse Engel, Kumar Krishna Agrawal, Shuo Chen, Ishaan Gulrajani, Chris Donahue, and Adam Roberts,  “Gansynth: Adversarial neural audio synthesis,”  arXiv preprint arXiv:1902.08710, 2019.
  * [16] Z. Kong, W. Ping, J. Huang, K. Zhao, and B. Catan,  “Diffwave: A versatile diffusion model for audio synthesis,”  arXiv preprint arXiv:2009.09761, 2020.
  * [17] Antoine Caillon and Philippe Esling,  “Rave: A variational autoencoder for fast and high-quality neural audio synthesis,”  arXiv preprint arXiv:2111.05011, 2021.
  * [18] Gwendal Le Vaillant, Thierry Dutoit, and Sébastien Dekeyser,  “Improving synthesizer programming from variational autoencoders latent space,”  in 2021 24th International Conference on Digital Audio Effects (DAFx). IEEE, 2021, pp. 276–283.
  * [19] Pritish Chandna, Merlijn Blaauw, Jordi Bonada, and Emilia Gómez,  “Wgansing: A multi-voice singing voice synthesizer based on the wasserstein-gan,”  in 2019 27th European signal processing conference (EUSIPCO). IEEE, 2019, pp. 1–5.
  * [20] Q. Huang, D. S. Park, T. Wang, T. I. Denk, A. Ly, N. Chen, Z. Zhang, Z. Zhang, J. Yu, C. Frank, et al.,  “Noise2music: Text-conditioned music generation with diffusion models,”  arXiv preprint arXiv:2302.03917, 2023.
  * [21] Siyuan Shan, Lamtharn Hantrakul, Jitong Chen, Matt Avent, and David Trevelyan,  “Differentiable wavetable synthesis,”  in ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022, pp. 4598–4602.
  * [22] Boris Kuznetsov, Julian D Parker, and Fabián Esqueda,  “Differentiable iir filters for machine learning applications,”  in Proc. Int. Conf. Digital Audio Effects (eDAFx-20), 2020, pp. 297–303.
  * [23] Marco A Martínez Ramírez, Oliver Wang, Paris Smaragdis, and Nicholas J Bryan,  “Differentiable signal processing with black-box audio effects,”  in ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2021, pp. 66–70.
  * [24] Christian J Steinmetz, Jordi Pons, Santiago Pascual, and Joan Serrà,  “Automatic multitrack mixing with a differentiable mixing console of neural audio effects,”  in ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2021, pp. 71–75.
  * [25] Christian J Steinmetz, Nicholas J Bryan, and Joshua D Reiss,  “Style transfer of audio effects with differentiable signal processing,”  arXiv preprint arXiv:2207.08759, 2022.
  * [26] Bo-Yu Chen, Wei-Han Hsu, Wei-Hsiang Liao, Marco A Martínez Ramírez, Yuki Mitsufuji, and Yi-Hsuan Yang,  “Automatic dj transitions with differentiable audio effects and generative adversarial networks,”  in ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022, pp. 466–470.
  * [27] N. Masuda and D. Saito,  “Improving semi-supervised differentiable synthesizer sound matching for practical applications,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, vol. 31, pp. 863–875, 2023.
  * [28] J. M. Chowning,  “The synthesis of complex audio spectra by means of frequency modulation,”  Journal of the Audio Engineering Society, vol. 21, no. 7, pp. 526–534, 1973.
  * [29] M. Lavengood,  “What makes it sound ’80s? the yamaha dx7 electric piano sound,”  Journal of Popular Music Studies, vol. 31, no. 3, pp. 73–94, 2019\.
  * [30] Eric Jang, Shixiang Gu, and Ben Poole,  “Categorical reparameterization with gumbel-softmax,”  arXiv preprint arXiv:1611.01144, 2016.


# Valence can control the nonexponential viscoelastic relaxation of
multivalent reversible gels

Hugo Le Roy †<EMAIL_ADDRESS>Université Paris-Saclay, CNRS, LPTMS, 91405,
Orsay, France  Institute of Physics, École Polytechnique Fédérale de Lausanne
(EPFL), 1015 Lausanne, Switzerland  Jake Song  These authors contributed
equally to this work Department of Materials Science and Engineering,
Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA
02139, USA  David Lundberg  Department of Chemical Engineering, Massachusetts
Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA
Aleksandr V. Zhukhovitskiy  Department of Chemistry, Massachusetts Institute
of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA  Department
of Chemistry, University of North Carolina at Chapel Hill; Chapel Hill, NC
27599, USA  Jeremiah A. Johnson  Department of Chemistry, Massachusetts
Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA
Gareth H. McKinley  Department of Mechanical Engineering, Massachusetts
Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA
Niels Holten-Andersen  Department of Materials Science and Engineering,
Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA
02139, USA  Martin Lenz<EMAIL_ADDRESS>Université
Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France  PMMH, CNRS, ESPCI Paris, PSL
University, Sorbonne Université, Université de Paris, F-75005, Paris, France

###### Abstract

Gels made of telechelic polymers connected by reversible crosslinkers are a
versatile design platform for biocompatible viscoelastic materials. Their
linear response to a step strain displays a fast, near-exponential relaxation
when using low valence crosslinkers, while larger supramolecular crosslinkers
bring about much slower dynamics involving a wide distribution of time scales
whose physical origin is still debated. Here, we propose a model where the
relaxation of polymer gels in the dilute regime originates from elementary
events in which the bonds connecting two neighboring crosslinkers all
disconnect. Larger crosslinkers allow for a greater average number of bonds
connecting them, but also generate more heterogeneity. We characterize the
resulting distribution of relaxation time scales analytically, and accurately
reproduce stress relaxation measurements on metal-coordinated hydrogels with a
variety of crosslinker sizes including ions, metal-organic cages, and
nanoparticles. Our approach is simple enough to be extended to any crosslinker
size and could thus be harnessed for the rational design of complex
viscoelastic materials.

## Introduction

Soft hydrogels are ubiquitous in biology and dictate the mechanics of cells
and tissues [1]. Due to their biocompatibility, synthetic hydrogels are thus
prime candidates to serve as robust soft tissue implants, although fine
control of their viscoelastic properties is crucial for their success in this
role [2, 3]. In simple viscoelastic materials, stress relaxes according to a
single exponential with a single relaxation time. This is not however the case
for most biological materials such as cells [4], tissues [5], mucus [6], and
biofilms [7]. Instead, their relaxation is characterized by a broad
distribution of relaxation times [8]. Such relaxation is often heuristically
described by a stretched exponential:

$\sigma(t)\propto e^{-(t/\tau)^{\alpha}},$ (1)

where smaller values of the exponent $\alpha\in]0;1[$ denote broader
distributions of relaxation time scales [9]. Other similarly phenomenological
fitting functions include power-law dependences of $\sigma$ on $t$ [10, 11,
12] and log-normal distributions of the relaxation times [13, 14].

Figure 1:  High-valence crosslinkers yield a slow, potentially complex
unbinding dynamics (a) Hydrogels held together by small crosslinkers relax
over the time scale associated with the unbinding of a single polymer strand.
There, stress decays exponentially in response to a step strain:
$\sigma(t)/\sigma(0)=\exp(-t/\tau)$. Since $t/\tau=\exp(\ln t-\ln\tau)$,
changes in the time scale $\tau$ shift the $\sigma$ vs. $\ln t$ relaxation
curve horizontally, but do not alter its shape. (b) In contrast, relaxation
events in the presence of high-valence crosslinkers require the simultaneous
unbinding of many polymer stands. The associated time scale is long and highly
variable depending on the number of strands involved in the “superbond” (grey
shade). As a result the stress relaxation of such gels is no longer
exponential and the precise shape of the relaxation curve strongly depends on
the valence of the crosslinkers.

Associative gels, which relax by a succession of binding and re-binding events
[15], offer a promising route to design controllable viscoelastic materials.
It is thus possible to tune their relaxation time by adjusting the chemical
binding energy of their crosslinkers [16, 17]. Although this chemistry-based
approach allows tuning of the overall stress relaxation time as illustrated in
Fig. 1(a), less is known about the different approaches to tune the shape of
the stress relaxation curve of reversible hydrogels. Accordingly, most
existing models for the relaxation of multivalent gels focus on regimes
dominated by a single relaxation time scale [18], leading to exponential
relaxation [19]. Control over the distribution of relaxation time scales could
however be achieved in synthetic hydrogels connected with multivalent
crosslinkers such as nanoparticles [20], metal-organic cages [21], clay [22],
and latex beads [23], which are known to exhibit nonexponential viscoelastic
relaxation. Here, we aim to elucidate this emergence of a wide distribution of
time scales in materials with high valence crosslinkers to enable the rational
design of complex gels [Fig. 1(b)]. Here we use the term “valence” to
designate the number of polymer strands that a crosslinker can bind, a
property sometimes also sometimes referred to as their “functionality” [24].

We propose that the emergence of a broad distribution of relaxation time
scales arises from microscopic events consisting of the severing of the
physical connection between two crosslinkers. We first propose a model where
this connection, hereafter termed “superbond”, breaks if all its constitutive
crosslinkers are detached at the same time [Fig. 1(b)]. We show that the
breaking time of a superbond increases exponentially with the number of
strands involved, consistent with previous observation [25]. As a result of
this strong dependence, small spatial heterogeneities in the polymer
concentration may result in widely different relaxation times from one
superbond to the next. Such exponential amplification of relaxation times
originating from small structural differences forms the basis of models
previously used to describe the relaxation of soft glasses [26, 27, 28]. In
contrast with these studies, our approach explicitly models the microscopic
basis of this amplification. That allows it to not only recover relaxation
curves virtually indistinguishable from those discussed in previous studies,
but to also predict the influence of temperature and crosslinker valence on
the macroscopic stress relaxation observed in the resulting gel. The details
of the polymer strand morphology are not central to this influence, and we
thus enclose them in a few effective parameters which could be derived from
first principles in specialized models related to specific implementations of
our basic mechanism. To confirm these predictions, we conduct experiments on
hydrogels with four distinct crosslinker types of different sizes, and find
that our model quantitatively reproduces multiple relaxation curves using this
small set of microscopic parameters. Finally, we show that several
phenomenological fitting functions used in the literature can be recovered as
asymptotic regimes of our analytical model.

## Results

### Model of a single superbond

We first model a single superbond in the simple, experimentally relevant [21]
case of strong interactions combined with short polymers, which implies
negligible entanglements. In the limit of very large and rigid crosslinkers,
the polymer layer around a crosslinker is locally planar, its structure is not
affected by small fluctuations of polymer concentration and its thickness
fixes the distance between crosslinkers [29, 21]. We thus use the simplifying
assumption that all individual bonds participating in a superbond are
identical and non-interacting, an approximation whose validity we ultimately
assess through comparisons with experiments.

We model the attachment and detachment of a single polymer strand from a pair
of crosslinkers as shown in Fig. 2(a). When both its ends are bound, the
strand may or may not connect two different crosslinkers. The corresponding
“bridging” and “looping” states have the same energy since we assume the
polymer strand to be completely flexible, and we denote by $\Delta S$ the
entropy difference between them. To transition between these two states, the
strand must disconnect one of its ends and form a “dangling” state. The
disconnection of the strand in this state implies an energy barrier $\Delta E$
that is much larger than the thermal energy $k_{B}T=\beta^{-1}$. This implies
that the dangling state is short-lived, and thus need not be explicitly
included in our modeling. Our approximation scheme implies that the overall
rate $\omega^{+}$ and $\omega^{-}$ to go from the looping to the bridging
state and back are constant. They read

$\omega^{+}=\frac{1}{\tau_{0}}e^{-\beta\Delta
E}\qquad\omega^{-}=\frac{1}{\tau_{0}}e^{-\beta\Delta E+\Delta S},$ (2)

where the typical time scale $\tau_{0}$ takes into account the entropy
difference between the looping and dangling state. At equilibrium, we denote
the probability for a single polymer strand to create a bridge as
$p_{\text{on}}=1-p_{\text{off}}=1/(1+e^{\Delta S})$.

Figure 2: We model superbond breaking as the disconnection of many independent
polymer strands. (a) Disconnecting a single polymer strand requires going
through a high-energy, short-lived dangling state (larger arrows indicate
faster transitions). The looping and bridging states both have two polymer-
crosslinker bonds, and therefore have the same energy. (b) Individual strands
in a superbond attach and detach independently, resulting in a one-dimensional
random walk in the number $n$ of attached strands [Eq. (3)]. Here we only draw
the bridging strands, not the looping strands.

We now consider the dynamics of a superbond involving $N$ polymer strands.
Within our approximation of independent attachment and detachment of the
individual polymer strands, the superbond undergoes the Markov process
illustrated in Fig. 2(b) and the probability $P_{n}(t)$ for $n$ strands to
create bridges between the two crosslinkers at time $t$ satisfies the master
equation

$\displaystyle\partial_{t}P_{n}(t)=$
$\displaystyle\,(N-n+1)\,\omega^{+}P_{n-1}(t)+(n+1)\,\omega^{-}P_{n+1}(t)$
$\displaystyle-[(N-n)\,\omega^{+}+n\,\omega^{-}]P_{n}(t),$ (3)

which ensures that the number of bridging polymers can never be greater than
$N$.

To determine the rate at which a superbond breaks, we set an absorbing
boundary condition $P_{0}(t)=0$ and define its survival probability as
$S(t)=\sum_{n=1}^{N}P_{n}(t)$. In the limit $N\gg 1$ where a large number of
strands are involved in the superbond, we show in Sec. S1 of the Supplementary
Information that the detachment of the two beads is analogous to a Kramers
escape problem. We thus prove that the survival probability decays as a single
exponential $S(t)=\exp(-t/\tau_{N})$ [30] with an average detachment time

$\tau_{N}\underset{N\rightarrow\infty}{\sim}\frac{\tau_{0}e^{\beta\Delta
E}}{Np_{\text{off}}^{N}}.$ (4)

The breaking of the superbond can thus be assimilated to a Poisson process
with rate $1/\tau_{N}$ regardless of the initial condition $P_{n}(0)$. The
strong, exponential dependence of $\tau_{N}$ on $N$ implies that any
dispersity in the number of strands involved in a superbond may result in a
wide distribution of time scales.

### Model of the relaxation of a gel

Two factors influence the dispersity of $N$. First, its value is constrained
by the available space at the surface of each crosslinker, which we model by
setting an upper bound $N_{\text{sat}}$ on the number of polymer strands (in a
loop or a bridge) participating in any superbond. Second, depending on the
local density of polymer in the vicinity of the superbond, the actual number
of strands present may be lower than $N_{\text{sat}}$. In the regime where the
polymer solution surrounding the crosslinkers is dilute, polymer strands are
independently distributed throughout the system. As a result, the distribution
of local strand concentrations within a small volume surrounding a superbond
follows a Poisson distribution. We thus assume that $N$ is also described by a
Poisson distribution up to its saturation at $N_{\text{sat}}$:

$p(N)=\begin{cases}\frac{\bar{N}^{N}e^{-\bar{N}}}{N!}&\text{for
}N<N_{\text{sat}}\\\
\sum_{K=N_{\text{sat}}}^{+\infty}\frac{\bar{N}^{K}e^{-\bar{N}}}{K!}&\text{for
}N=N_{\text{sat}}\end{cases},$ (5)

where $\bar{N}$ would be the average number of strands in a superbond in the
absence of saturation and thus depends on the ratio of polymer to crosslinker
concentration. Note that the specific form of the distribution used in Eq. (5)
does not significantly modify our results, as discussed later.

In response to a step strain, we assume that each superbond is stretched by an
equal amount and resists the deformation with an equal force prior to
breaking. Superbonds may subsequently reform, but the newly formed bonds are
not preferentially stretched in the direction of the step strain and therefore
do not contribute to the macroscopic stress on average. Denoting by $t=0$ the
time at which the step strain is applied and by $\sigma(t)$ the resulting
time-dependent shear stress, the progressive breaking of the initial
superbonds results in the following stress response function:

$\frac{\sigma(t)}{\sigma(t=0)}=\sum_{N=1}^{N_{\text{sat}}}\frac{p(N)}{1-p(0)}e^{-t/\tau_{N}}.$
(6)

While the breaking times $\tau_{N}$ are unaffected by the applied stress in
the linear response regime, nonlinearities could easily be included in our
formalism by making $\Delta S$ stress-dependent and thus favor strand
detachment. The relaxation described in Eq. (6) occurs in two stages. At long
times $t\gg\tau_{N_{\text{sat}}}$, few short-lived superbonds remain.
Saturated superbonds ($N=N_{\text{sat}}$) dominate the response and Eq. (6) is
dominated by the last term of its sum. As a result, the stress relaxes
exponentially over time, as seen from the linearity of the log-lin curves of
Fig. 3(a) for large values of $t$. Systems with smaller values of
$N_{\text{sat}}$ manifest this regime at earlier times; in the most extreme
case, the relaxation of a system where superbonds involve at most a single
polymer strand ($N_{\text{sat}}=1$) is fully exponential and extremely fast as
compared to systems with higher $N_{\text{sat}}$. Over short times
($t\ll\tau_{N_{\text{sat}}}$), stress relaxation involves multiple time
scales. This nonexponential regime is apparent on the left of Fig. 3(a). These
two regimes have already been reported in several experimental gels connected
by multivalent crosslinkers [1, 29].

Figure 3: (a) Disperse, high-valence superbonds initially display a
nonexponential mechanical relaxation, then cross over to an exponential regime
when only the saturated superbonds remain. Curves plotted from Eq. (6) with
$p_{\text{off}}=0.2$, $\bar{N}=10$ and different values of $N_{\text{sat}}$ as
indicated on each curve. (b) Relationship between the stretch exponent
$\alpha$ quantifying the nonexponential character of the relaxation and the
microscopic parameter $N_{\text{sat}}/\bar{N}$. Here $p_{\text{off}}=0.2$. A
low $N_{\text{sat}}/\bar{N}$ gives an exponential relaxation ($\alpha\simeq
1$), while a larger $N_{\text{sat}}/\bar{N}$ leads to a more complex behavior
($\alpha<1$). While $\alpha$ appears to converge to a finite value for large
$N_{\text{sat}}/\bar{N}$ for the largest values of $\bar{N}$, this behavior is
contingent on our choice of fitting interval. This issue does not affect the
rest of the curves. Large stars correspond to the curves represented in (a).
Inset: illustration of the quality of the fits between the heuristic stretched
exponential [dashed orange line, Eq. (1)] and our prediction [solid blue line,
Eq. (6)]. Figure 4: Stress relaxation function for four experimental systems
with increasing crosslinker valences (see Table 1 for values). Here we use a
log-lin scale [unlike in Fig. 3(a)] to facilitate the visualization a large
range of time scales. Alternate representations are available as Fig. S5.
Symbols are experimental datapoints, and the lines are the associated fitting
curves. Insets: time-temperature collapsed data obtained by a rescaling
$t\rightarrow te^{\beta\Delta E}$.

While Eq. (6) is not identical to the stretched exponential of Eq. (1), the
inset of Fig. 3(b) shows that they are remarkably close in practice. We thus
relate the stretch exponent $\alpha$ to the saturation number $N_{\text{sat}}$
by fitting a stretched exponential to our predicted stress response function
over the time interval required to relax $90\%$ of the initial stress [Fig.
3(b)]. The fits are very close matches, and consistently give correlation
factors $r^{2}>0.98$ (see detailed plots in Supplementary Fig. S2). If
$N_{\text{sat}}\lesssim 0.5\bar{N}$ then $\alpha\simeq 1$, indicating a
nearly-exponential relaxation. Indeed, in that case superbond saturation
occurs well before the peak of the Poisson distribution of $N$. Physically,
this implies that the local polymer concentration surrounding most superbonds
is sufficient to saturate them. As almost all superbonds are saturated, they
decay over the same time scale $\tau_{N_{\text{sat}}}$. As a result, the
material as a whole displays an exponential relaxation. For larger values of
$N_{\text{sat}}$, the Poisson distribution is less affected by the saturation
and the dynamics is set by the successive decay of superbonds involving an
increasing number of strands, implying lower values of $\alpha$. The larger
the value of $\bar{N}$, the sharper the crossover between these two regimes.

### Experiments

To validate our model of the effect of crosslinker valence on hydrogel
relaxation, we perform step-strain experiments of poly(ethylene glycol)
(PEG)-based gels involving a range of crosslink valences and compare Eq. (6)
to the resulting relaxation curves. We implement strand-crosslink bonds based
on two sets of metal-cooordination chemistry: nitrocatechol-$\text{Fe}^{3+}$
coordination and bispyridine-$\text{Pd}^{2+}$ coordination. The first set of
gels is made with nitrocatechol-functionalized PEG crosslinked by single
$\text{Fe}^{3+}$ ions with an estimated valence of 3, and by iron oxide
nanoparticles. The nanoparticles have a mean diameter of $7\,\text{nm}$ with a
surface area that allows a valence of $\sim 100$ ligands [20]. The second set
of gels are made with bispyridine-functionalized PEG, wherein
bis-$meta$-pyridine ligands induce self-assembly of gels that are crosslinked
by $\text{Pd}_{2}$$\text{L}_{4}$ nanocages with a valence of 4, and
bis-$para$-pyridine ligands induce self-assembly of gels that are crosslinked
by $\text{Pd}_{12}$$\text{L}_{24}$ nanocages with a valence of 24 [31]. As
shown in Fig. 4, these four distinct gel designs result in a broad range of
relaxation behaviors. Overall, large-valence gels and lower temperatures
result in longer relaxation times, consistent with the illustration of Fig.
1(a). The relaxation curves associated to high-valence crosslinkers are also
less steep, consistent with the involvement of a broader distributions of
relaxation times and the schematic of Fig. 1(b).

crosslinker | $\text{Fe}^{3+}$ | Pd2L4 | Pd12L24 | nanoparticles
---|---|---|---|---
estimated valence | $3$ | $4$ | $24$ | $100$
$\bm{N_{\text{sat}}}$ | $1$ | $4$ | $7$ | $17$
$\bm{\Delta E}$ (units of $k_{B}T$) | $28$ | $24$ | $24$ | $24$
$\bm{p}_{\mathbf{off}}$ | $0.08$ | $0.15$ | $0.15$ | $0.36$
$\bm{\tau_{1}}\text{ at $T=300\,\text{K}$}$ (s) | $1.7$ | $6.1$ | $1.9$ | $0.1$
$\bm{\bar{N}}$ | $1$ | $6$ | $9$ | $14$
Table 1: Estimated and fitted parameters involved in the comparison between
experiment and theory in Fig. 4. The energies are given in units of $k_{B}T$
for $T=300\,\text{K}$. Instead of displaying the parameter $\tau_{0}$, we
present the more easily interpreted unbinding time of a single polymer strand
at $300K$, namely $\tau_{1}=\tau_{0}e^{\beta_{300}\Delta E}/p_{\text{off}}$.

At a more detailed level, our assumption that the dynamics of single polymer
strand proceeds independently of its environment implies the existence of a
single energy scale $\Delta E$. As a result, we predict that all time scales
involved in the relaxation are proportional to $\exp(-\beta\Delta E)$. We
confirm this through a time-temperature collapse shown in the insets of Fig. 4
(see Sec. S4 of the Supplementary Information for details). This collapse
provides us with the value of $\Delta E$ for each of our four systems, which
we report in Table 1. The binding energy value $\Delta E$ for the
$\text{Pd}_{2}$$\text{L}_{4}$ and $\text{Pd}_{12}$$\text{L}_{24}$ gels match,
as expected from the fact that they originate from the same composition. On
the other hand, the $\Delta E$ of the nanoparticles gel cannot be directly
compared to the $\text{Fe}^{3+}$ ion gel because of their dramatic physico-
chemical properties as well as the acidic pH used to facilitate mono-
coordination in nanoparticles.

To compare the temperature-collapsed curves to our prediction of Eq. (6), we
fit the parameters $p_{\text{off}}$, $\tau_{0}$, $\bar{N}$ and
$N_{\text{sat}}$ across multiple temperatures. The resulting fits, shown in
Fig. 4, display a good agreement between the theory and experiments across up
to 4 orders of magnitude in time scales. The fitted values of $N_{\text{sat}}$
are moreover consistent with the estimated valence of the crosslinkers (Table
1), which confirms our interpretation of the physical origin of
$N_{\text{sat}}$. A more direct relation between the valence and
$N_{\text{sat}}$ would require additional information about the structure of
each category of gels, and in particular the number of nearest neighbors of an
individual crosslink, which remains debated for these types of gels [32, 33,
34]. The possible clustering of crosslinkers may also influence this relation,
as nanocage systems similar to ours [35] may form higher-order nanocage
structures. While such complexities as well as possible imperfections in our
fitting procedure complicate the literal interpretation of the fitted values
of $N_{\text{sat}}$, these nonetheless clearly increase with increasing
crosslinker valence. The fit also supports the notion that the mean number of
strands per superbond $\bar{N}$ accounts for the distribution of relaxation
time scales in our gels. The $\text{Fe}^{3+}$ gel thus displays an exponential
relaxation consistent with $\bar{N}=N_{\text{sat}}=1$. The higher-valence
Pd2L4 and Pd12L24 systems have a complex relaxation at early times followed by
an exponential behavior, as expected for $\bar{N}\simeq N_{\text{sat}}>1$. As
expected from our model, the crossover time $\tau_{N_{\text{sat}}}$ separating
the two regimes is larger in the higher-valence Pd12L24 gel. Finally, the
high-valence nanoparticle system shows an extended complex relaxation
associated with $\bar{N}<N_{\text{sat}}$, thus confirming that all the
qualitative relaxation regimes discussed in the previous sections are
experimentally relevant.

### Distribution of relaxation time scales

To further visualize the differences between the responses of our gels, we
plot the distributions of relaxation times $p(\tau)$ for our fitted model in
Fig. 5. The Fe3+ gels, which relax according to a single exponential and whose
$p(\tau)$ are therefore delta functions, are not represented there. In the
Pd2L4 and Pd12L24 systems, a distribution characterized by an initially
decreasing distribution of time scales is interrupted by a valence-dependent
maximum relaxation time $\tau_{N_{\text{sat}}}$. That time is comprised within
the range of time scales observed in Fig. 4, accounting for the crossover to
an exponential relaxation within this range. In nanoparticle systems, by
contrast, the crossover occurs much later and thus cannot be directly observed
in experiments. In all cases, the precise form of the distribution of time
scales used in the domain $\tau<\tau_{N_{\text{sat}}}$ does not critically
affect the predicted relaxation curves. Indeed, we show in Sec. S6 of the
Supplementary Information that replacing the Poisson distribution of Eq. (5)
with other distributions with the same mean and variance lead to essentially
indistinguishable predictions over experimentally observable time scales. This
emphasizes the robustness of our predictions to the details of that choice of
distribution. They are instead primarily determined by the mean and maximum
superbond sizes, $\bar{N}$ and $N_{\text{sat}}$.

Figure 5: Distribution of relaxation times corresponding to the theoretical
plots of Fig. 4. The distribution associated with the Fe3+ gel is a delta
function for $\tau/\tau_{1}=1$ and is thus not represented on this graph. The
time distributions are given by a Poisson distribution cut for
$\tau=\tau_{N_{\text{sat}}}$ (vertical lines), as described in Eq. (6). The
dotted lines represent what these distributions would be in the absence of
this saturation.

In the limit of large $\bar{N}$ and even larger $N_{\text{sat}}$, the complex
relaxation phase of our model characterized by $t<\tau_{N_{\text{sat}}}$ may
display analytical behaviors identical to some widely used rheological fitting
functions. In this regime, the Poisson distribution $p(N)$ of Eq. (5) goes to
a Gaussian. Since according to Eq. (4) the variable $N$ is essentially the
logarithm of the relaxation time $\tau$ for $\bar{N}\gg 1$, this results in a
log-normal distribution of relaxation time scales:

$\displaystyle p(\tau)=$
$\displaystyle\frac{1}{\sqrt{2\pi\bar{N}}|\ln{p_{\text{off}}|}\tau}$ (7)
$\displaystyle\times\exp\left\\{-\frac{\left[\ln\tau+\ln{\frac{\bar{N}}{\tau_{1}}}+\ln
p_{\text{off}}\left(\bar{N}-1\right)\right]^{2}}{2\bar{N}(\ln{p_{\text{off}}})^{2}}\right\\}.$

This result adds additional insights to this widely used fitting functional
form, as it allows to relate the mean and variance of the distribution to the
underlying crosslinker-scale parameters [13, 14]. It moreover offers a
potential molecular-level justification for its use in describing the complex
relaxation of systems with multivalent crosslinks. In the alternative case
where $p(N)$ is a decaying exponential, our model results in power-law
distributed relaxation time scales and the stress response function takes the
form

$\sigma(t)\propto t^{-\gamma},\quad\text{with}\quad\gamma=\frac{1}{\bar{N}|\ln
p_{\text{off}}|}.$ (8)

This result may also be presented in terms of the dependence of the storage
and loss moduli on the frequency $\omega$ in an oscillatory rheology
experiment. We thus predict that for $\gamma<1$

$\displaystyle G^{\prime}(\omega)$
$\displaystyle\approx\begin{cases}\omega^{2}&\text{for
}\omega\ll\tau_{N_{\text{sat}}}^{-1}\\\ \omega^{\gamma}&\text{for
}\tau_{N_{\text{sat}}}^{-1}\ll\omega\ll\tau_{1}^{-1}\\\ \omega^{0}&\text{for
}\tau_{1}^{-1}\ll\omega\end{cases}$ (9a) $\displaystyle
G^{\prime\prime}(\omega)$
$\displaystyle\approx\begin{cases}\omega^{1}&\text{for
}\omega\ll\tau_{N_{\text{sat}}}^{-1}\\\ \omega^{\gamma}&\text{for
}\tau_{N_{\text{sat}}}^{-1}\ll\omega\ll\tau_{1}^{-1}\\\ \omega^{-1}&\text{for
}\tau_{1}^{-1}\ll\omega\end{cases}$ (9b)

The results larger values of $\gamma$ and detailed derivations of Eqs. (7-9)
are shown in the Supplementary Information. Again, this result has the
potential to account for the power-law relaxation observed in many rheological
systems [10, 11, 12], in addition to providing a link to their microscopic
constituents. Overall, these results suggest a possible control of the
system’s rheology through the characteristics of $p(N)$, which could in turn
be modulated through the spatial distribution of the polymer strands and the
dispersity of the crosslinkers.

## Discussion

Our simple model recapitulates a wide range of rheological behaviors in
multivalent systems based on two key superbond parameters: the mean size
$\bar{N}$ and the maximum size $N_{\text{sat}}$. These respectively control
the amplitude of the fluctuations in superbond size and the longest superbond
relaxation time scale. Prior to the longest relaxation time, the system
displays an increasingly nonexponential response for increasing $\bar{N}$
[Fig. 3(b)]. Beyond it, it crosses over into exponential relaxation. In
contrast with widely used phenomenological fitting parameters, our two
variables yield reliable insights into the underlying microscopic dynamics, as
demonstrated by the agreement of their fitting values with our a priori
knownledge of four experimental systems covering a wide range of values of
$\bar{N}$ and $N_{\text{sat}}$.

Our model bears a mathematical similarity with standard random energy trap
models [36]. There, a long-tailed relaxation emerges from a short-tailed
distribution of trap depths due to the exponential dependence of the
relaxation times on the trap depths. Similarly, here a nonexponential
relaxation emerges from a short-tailed distribution of superbond sizes $N$
[Eq. (5)] thanks to the exponential dependence of $\tau_{N}$ on $N$ [Eq. (4)].
In contrast with trap models, however, our model does not predict a glass
transition upon a lowering of temperature. It instead displays a simple
Arrhenius time-temperature relation, consistent with the experimental
collapses in the insets of Fig. (4). Here again, an additional benefit of our
approach is the direct connection between the predicted relaxation and
experimentally accessible parameters such as the crosslinker surface (through
$N_{\text{sat}}$).

Our model’s focus on the collective aspects of superbond breaking and the
characteristics of the crosslinkers implies that it encloses most of the
physics of the polymer strands within a few mesoscopic parameters, mainly
$\tau_{0}$ and $\Delta S$. Within our approach, the morphology of the polymer
thus does not affect the form of our relaxation, although it may lead to a
rescaling of the relaxation times of Eq. (4). This formulation remains valid
as long as the length and concentration of the polymer strands is low enough
that the polymer strands do not become significantly entangled, which could
spoil the Poissonian attachment/detachment process of Eq. (2). Even in this
case however this equation may not be strictly valid, as the polymer layer in
a superbond with many bound crosslinkers tends to be more compressed than in
one with few. This effect should leads to a smooth (likely power law)
dependence of $\omega^{\pm}$ on $N$, which would preserve the dominance of the
much more abrupt exponential dependence of $\tau_{N}$ on $N$. As a result,
while such polymer brush effects could induce corrections in our estimations
of the model parameters, the basic mechanism outlined here should still hold
in their presence.

Our model reproduces several qualitative characteristics of the rheology of
multivalent gels, such as the strong influence of the crosslinker valence,
Arrhenius temperature dependence and the transition between a nonexponential
and an exponential regime at long times. Due to its simple, widely applicable
microscopic assumptions, we believe that it could help shed light on and
assist the design of a wide range of multivalent systems. Beyond composite
gels, it could thus apply to RNA-protein biocondensates where multivalent
interactions between proteins are mediated by RNA strands [37], as well as
cytoskeletal systems where filaments linked to many other filaments display a
slow relaxation reminiscent of that of our multivalent crosslinkers [38].

## Methods

### Materials

4-arm poly(ethylene glycol) bis(acetic acid N-succinimidyl) ester (4-arm PEG-
NHS) (MW = 10 kDa), 1-arm poly(ethylene glycol) bis(acetic acid
N-succinimidyl) ester (1-arm PEG-NHS) (MW = 2000 Da) were purchased from
JenKem Technology. Sodium sulfate ($\text{Na}_{2}\text{SO}_{4}$), sodium
nitrite ($\text{NaNO}_{2}$), iron(III) acetylacetonate
$(\text{Fe(acac)}_{3})$, hydrochloric acid (HCl), dopamine hydrochloride,
triethylamine (TEA), N-methylmorpholine (NMM), dimethyl sulfoxide (DMSO),
methanol (MeOH), ethanol (EtOH), dichloromethane (DCM), N, N-dimethylformamide
(DMF), diethyl ether ($\text{Et}_{2}\text{O}$), and chloroform
($\text{CHCl}_{3}$), were purchased from Sigma-Aldrich. All chemicals were
used without further purification.

### Synthesis of 1-arm PEG-Catechol (C)

228 mg dopamine hydrochloride is neutralized for 15 min with 0.3 mL NMM in 7.5
mL dry DMF under $\text{N}_{2}$ atmosphere. Then 1 g mPEG-NHS (MW = 2000 Da)
dissolved in 7.5 mL DMF is added, and the mixture is stirred with N2
protection at room temperature for 24 hours. The reacted solution is acidified
by adding 15 mL 1M HCl (aq), and the product is extracted with
$\text{CHCl}_{3}$ 3 times. The organic layers are pooled together, dried with
$\text{NaSO}_{4}$, and solvent is removed by rotary evaporation. Finally, the
product concentrate is precipitated in cold $\text{Et}_{2}\text{O}$ (-20 ∘C),
filtered and dried. 1H NMR (300 MHz, D2O) $\delta$ (ppm): 6.7-6.8 (m, 3H,
aromatic), 3.3-4.0 (m, $\text{-O-CH}_{2}\text{-CH}_{2}\text{-}$), 3.4 (t, 2H,
CH2 adjacent to aromatic ring), 2.7 (t, 2H, -CH2-NH-CO-).

### Synthesis of 4-arm PEG-Nitrocatechol (NC)

178 mg of nitrodopamine hydrogen sulfate is neutralized for 15 min with 110
$\mu$L NMM in 4 mL dry DMF under N2 atmosphere. Then 1 g 4-arm PEG-NHS (MW =
10 kDa) dissolved in 4 mL DMF is added, and the mixture is stirred with
$\text{N}_{2}$ protection at room temperature for 24 hours. The reacted
mixture is mixed with 15 mL 1M HCl(aq), dialyzed with water (MWCO = 3500 Da)
for 2 days (water exchanged for more than 5 times), and freeze-dried. 1H NMR
(300 MHz, D2O) $\delta$ (ppm): 7.6 (m, 1H, aromatic), 6.7 (m, 1H, aromatic),
3.6-3.9 (m, -O-CH2-CH2-), 3.5 (t, 2H, CH2 adjacent to aromatic ring), 3.1 (t,
2H, -CH2-NH-CO-).

### Synthesis of Fe3O4 nanoparticles (NPs)

Bare Fe3O4 NPs are synthesized following previously reported methods (Li et
al. 2016). 100 mg as-synthesized NPs are re-dispersed in 80 mL of 1:1 (v/v)
solution of CHCl3 and DMF, and 100 mg 1-arm PEG-C is added. The mixture is
homogenized and equilibrated by pulsed sonication (pulse: 10 s on + 4 s off;
power: 125 W) for 1 hour. Then the mixture is centrifuged at 10000 rpm for 10
min to remove any aggregates, and rotary evaporated at 50 ∘C, 30 mbar to
remove CHCl3. Then the NP solution is precipitated in 150 mL cold Et2O (-20
∘C). The precipitate is re-dispersed in H2O, and freeze-dried. The resulting
NPs are 7 nm in diameter.

### Preparation of the Fe$3^{+}$-NC gels

Preparation procedure is similar to a previously reported protocol [39],
except that the gel is made in DMSO instead of H2O. $50\,\mu$L of 200 mg/mL
4-arm PEG-NC solution in DMSO is mixed with $16.7\,\mu$L of 80 mM FeCl3
solution in DMSO (ligand: Fe3+ molar ratio of 3:1). Then $33.3\,\mu$L DMSO and
$13.8\,\mu$L TEA is added to facilitate deprotonation, and a gel is formed.

### Preparation of the Pd2L4 gels

The synthesis of polymer and gel preparation procedures for P2L4 is the same
as a reported protocol [31] with minor modifications. The annealing of the
Pd2L4 polyMOC gel was done at 60°C for 1 hr instead of 80°C for 4 hr, and 1.05
equivalent of Pd(NO3)2. 2 H2O (relative to bifunctional polymer ligand) was
used instead of 1 equivalent.

### Preparation of the Pd12L24 gels

The synthesis of polymer and gel preparation procedures for polyMOC is the
same as a reported protocol [31]

### Preparation of the NP gels

Preparation procedure is the same as the reported protocol [40]. Briefly,
$20\,\text{mg}$ PEGylated Fe3O4 NPs (equivalent to $20\,\text{mg}$ Fe3O4 core)
and $20\,\text{mg}$ 4-arm PEG-NC are mixed in a $0.2\,M$ HCl aqueous solution.
The solution mixture ($\text{pH}=2$) is transferred into a mold and sealed,
and a solid gel is obtained after curing in a 50 ∘C oven for 24 hours.

### Rheology

Stress relaxation measurements are done on an Anton Paar rheometer with
parallel plate geometry ($10\,\text{mm}$ diameter flat probe for NP gels and
polyMOC gels, and $25\,\text{mm}$ diameter cone probe for Fe3+ gels). All
tests are done immediately after transferring the gel sample onto the sample
stage. A Peltier hood is used for all experiments to control the measurement
temperature and prevent solvent evaporation. H2O based samples are furthermore
sealed with mineral oil before experimentation to reduce the evaporation rate.
Relaxation tests were performed by applying a $\gamma=0.005$. step strain for
the NP gel, and $\gamma=0.02$ step strain for the other three systems.

###### Acknowledgements.

HLR and ML thank Thibaut Divoux for valuable comments on the manuscript. This
work was supported by Marie Curie Integration Grant PCIG12-GA-2012-334053,
“Investissements d’Avenir” LabEx PALM (ANR-10-LABX-0039-PALM), ANR grants
ANR-15-CE13-0004-03, ANR-21-CE11-0004-02, ANR-22-ERCC-0004-01 and
ANR-22-CE30-0024-01, as well as ERC Starting Grant 677532 to ML. ML’s group
belongs to the CNRS consortium AQV.

## Author contributions

Hugo Le Roy and Martin Lenz designed the theoretical model and performed the
theoretical analysis. Hugo Le Roy analyzed the data. Jake Song synthetized the
ion and nanoparticle gels and carried out their rheological characterization
under the supervision of Niels Holten-Andersen and Gareth H. McKinley.
Aleksandr Zhukhovitskiy and David Lunderberg synthetized the nanocage gels and
carried out their rheological characterization under the supervision of
Jeremiah A. Johnson. The article was written by Hugo Le Roy, Jake Song and
Martin Lenz. All others read and edited the manuscript.

## References

  * Ronsin  _et al._ [2016] O. Ronsin, C. Caroli, and T. Baumberger, The Journal of Chemical Physics 144, 064904 (2016).
  * Borges  _et al._ [2011] A. C. Borges, C. Eyholzer, F. Duc, P.-E. Bourban, P. Tingaut, T. Zimmermann, D. P. Pioletti, and J.-A. E. Månson, Acta Biomaterialia 7, 3412 (2011).
  * Chaudhuri  _et al._ [2015] O. Chaudhuri, L. Gu, M. Darnell, D. Klumpers, S. A. Bencherif, J. C. Weaver, N. Huebsch, and D. J. Mooney, Nature Communications 6, 6365 (2015).
  * Kollmannsberger and Fabry [2011] P. Kollmannsberger and B. Fabry, Annual review of materials research 41, 75 (2011).
  * Chaudhuri  _et al._ [2020] O. Chaudhuri, J. Cooper-White, P. A. Janmey, D. J. Mooney, and V. B. Shenoy, Nature 584, 535 (2020).
  * Celli  _et al._ [2009] J. P. Celli, B. S. Turner, N. H. Afdhal, S. Keates, I. Ghiran, C. P. Kelly, R. H. Ewoldt, G. H. McKinley, P. So, S. Erramilli, _et al._ , Proceedings of the National Academy of Sciences 106, 14321 (2009).
  * Fabbri and Stoodley [2016] S. Fabbri and P. Stoodley, The perfect slime—microbial extracellular polymeric substances , 153 (2016).
  * Song  _et al._ [2022] J. Song, N. Holten-Andersen, and G. H. McKinley, arXiv preprint arXiv:2211.06438 (2022).
  * Bouchaud [2008] J.-P. Bouchaud, in  _Anomalous Transport_ (John Wiley & Sons, Ltd, 2008) Chap. 11, pp. 327–345.
  * Blair and Burnett [1959] G. S. Blair and J. Burnett, British Journal of Applied Physics 10, 15 (1959).
  * Keshavarz  _et al._ [2017] B. Keshavarz, T. Divoux, S. Manneville, and G. H. McKinley, ACS Macro Letters 6, 663 (2017).
  * Balland  _et al._ [2006] M. Balland, N. Desprat, D. Icard, S. Féréol, A. Asnacios, J. Browaeys, S. Hénon, and F. Gallet, Physical Review E 74, 021911 (2006).
  * Epstein  _et al._ [2019] E. S. Epstein, L. Martinetti, R. H. Kollarigowda, O. Carey-De La Torre, J. S. Moore, R. H. Ewoldt, and P. V. Braun, Journal of the American Chemical Society 141, 3597 (2019).
  * Masurel  _et al._ [2015] R. J. Masurel, S. Cantournet, A. Dequidt, D. R. Long, H. Montes, and F. Lequeux, Macromolecules 48, 6690 (2015).
  * Alberto Parada and Zhao [2018] G. Alberto Parada and X. Zhao, Soft Matter 14, 5186 (2018).
  * Grindy  _et al._ [2015] S. C. Grindy, R. Learsch, D. Mozhdehi, J. Cheng, D. G. Barrett, Z. Guan, P. B. Messersmith, and N. Holten-Andersen, Nature Materials 14, 1210 (2015).
  * Rosales and Anseth [2016] A. M. Rosales and K. S. Anseth, Nature Reviews Materials 1, 1 (2016).
  * Semenov  _et al._ [1995] A. N. Semenov, J.-F. Joanny, and A. R. Khokhlov, Macromolecules 28, 1066 (1995).
  * Michel  _et al._ [2001] E. Michel, J. Appell, F. Molino, J. Kieffer, and G. Porte, Journal of Rheology 45, 1465 (2001).
  * Li  _et al._ [2016] Q. Li, D. G. Barrett, P. B. Messersmith, and N. Holten-Andersen, ACS nano 10, 1317 (2016).
  * Zhukhovitskiy  _et al._ [2016a] A. V. Zhukhovitskiy, J. Zhao, M. Zhong, E. G. Keeler, E. A. Alt, P. Teichen, R. G. Griffin, M. J. A. Hore, A. P. Willard, and J. A. Johnson, Macromolecules 49, 6896 (2016a).
  * Wang  _et al._ [2010] Q. Wang, J. L. Mynar, M. Yoshida, E. Lee, M. Lee, K. Okuro, K. Kinbara, and T. Aida, Nature 463, 339 (2010).
  * Chatterjee  _et al._ [2014] T. Chatterjee, A. I. Nakatani, and A. K. Van Dyk, Macromolecules 47, 1155 (2014).
  * Flory [1953] P. J. Flory,  _Principles of polymer chemistry_ (Cornell university press, 1953).
  * Gomez-Casado  _et al._ [2011] A. Gomez-Casado, H. H. Dam, M. D. Yilmaz, D. Florea, P. Jonkheijm, and J. Huskens, Journal of the American Chemical Society 133, 10849 (2011).
  * Bertin and Bouchaud [2003] E. M. Bertin and J.-P. Bouchaud, Physical Review E 67, 026128 (2003).
  * Sollich  _et al._ [1997] P. Sollich, F. Lequeux, P. Hébraud, and M. E. Cates, Physical Review Letters 78, 2020 (1997).
  * Trachenko and Zaccone [2021] K. Trachenko and A. Zaccone, Journal of Physics: Condensed Matter 33, 315101 (2021).
  * Raspaud  _et al._ [1996] E. Raspaud, D. Lairez, M. Adam, and J.-P. Carton, Macromolecules 29, 1269 (1996).
  * Texier [2000] C. Texier, Journal of Physics A: Mathematical and General 33, 6095 (2000).
  * Zhukhovitskiy  _et al._ [2016b] A. V. Zhukhovitskiy, M. Zhong, E. G. Keeler, V. K. Michaelis, J. E. P. Sun, M. J. A. Hore, D. J. Pochan, R. G. Griffin, A. P. Willard, and J. A. Johnson, Nature Chemistry 8, 33 (2016b).
  * Hsiao  _et al._ [2012] L. C. Hsiao, R. S. Newman, S. C. Glotzer, and M. J. Solomon, Proceedings of the National Academy of Sciences 109, 16029 (2012).
  * Colombo and Gado [2014] J. Colombo and E. D. Gado, Soft Matter 10, 4003 (2014).
  * Gu  _et al._ [2018] Y. Gu, E. A. Alt, H. Wang, X. Li, A. P. Willard, and J. A. Johnson, Nature 560, 65 (2018).
  * Li  _et al._ [2008] D. Li, J. Zhang, K. Landskron, and T. Liu, Journal of the American Chemical Society 130, 4226 (2008).
  * Bouchaud [1992] J. P. Bouchaud, Journal de Physique I 2, 1705 (1992).
  * Choi  _et al._ [2020] J.-M. Choi, A. S. Holehouse, and R. V. Pappu, Annual Review of Biophysics 49, 107 (2020).
  * Lieleg  _et al._ [2011] O. Lieleg, J. Kayser, G. Brambilla, L. Cipelletti, and A. R. Bausch, Nature Materials 10, 236 (2011).
  * Holten-Andersen  _et al._ [2011] N. Holten-Andersen, M. J. Harrington, H. Birkedal, B. P. Lee, P. B. Messersmith, K. Y. C. Lee, and J. H. Waite, Proceedings of the National Academy of Sciences 108, 2651 (2011).
  * Song  _et al._ [2020] J. Song, M. H. Rizvi, B. B. Lynch, J. Ilavsky, D. Mankus, J. B. Tracy, G. H. McKinley, and N. Holten-Andersen, ACS Nano 14, 17018 (2020).