# Combining Stochastic Explainers and Subgraph Neural Networks can Increase
Expressivity and Interpretability

Indro Spinelli 1 and Michele Guerra 2 and Filippo Maria Bianchi 2,3 and Simone
Scardapane 4
1 - Sapienza University of Rome - DI
Via Salaria 113 This work was partially supported by the CHIST-ERA grant
CHIST-ERA-19-XAI-009 and PNRR MUR project PE0000013-FAIR. 00198 Rome - Italy

2 - UiT the Arctic University of Norway - Department of Mathematics and
Statistics
Hansine Hansens veg 18  9019 Tromsø - Norway

3 - NORCE Norwegian Research Centre
Nygårdstangen 22  5838 Bergen - Norway

4 - Sapienza University of Rome - DIET
Via Eudossiana 18  00184 Rome - Italy

###### Abstract

Subgraph-enhanced graph neural networks (SGNN) can increase the expressive
power of the standard message-passing framework. This model family represents
each graph as a collection of subgraphs, generally extracted by random
sampling or with hand-crafted heuristics. Our key observation is that by
selecting “meaningful” subgraphs, besides improving the expressivity of a GNN,
it is also possible to obtain interpretable results. For this purpose, we
introduce a novel framework that jointly predicts the class of the graph and a
set of explanatory sparse subgraphs, which can be analyzed to understand the
decision process of the classifier. We compare the performance of our
framework against standard subgraph extraction policies, like random node/edge
deletion strategies. The subgraphs produced by our framework allow to achieve
comparable performance in terms of accuracy, with the additional benefit of
providing explanations.

##  1 Introduction

Graph neural networks (GNNs) are neural network models designed to adapt and
perform inference on graph domains [1]. While a few models were already
proposed in-between 2005 and 2009 [2, 3], the interest in GNNs has increased
dramatically over the last few years, thanks to the broader availability of
data, processing power, and automatic differentiation frameworks. Now, GNNs
are the state-of-the-art solution in a comprehensive set of scenarios.
Nevertheless, regulations require high task performance and a transparent
decision process [4].

For this reason, several researchers have investigated techniques to explain
GNNs’ predictions, primarily identifying the most critical portions of the
graph that contributed to producing a particular inference. The vast majority
of these techniques [5, 6] provide a post-hoc explanation, thus inferring the
reasons that led to a specific outcome by a trained model. However, recent
efforts toward “explainable-by-design” GNNs rather than post-hoc explainers
are opening up new, interesting approaches. For example, in [7], the authors
introduce an explainability term to let the network converge to a
“interpretable” local minimum to facilitate the work of post-hoc explanation
algorithms. Solutions like [8, 9] completely discard the notion of post-hoc
algorithms providing explanations directly in the main model output.

On a separate line of research, recent studies [10, 11, 12] demonstrated that
by providing the GNN with subgraphs that give different views of the same
graph, it is possible to increase the expressive power of the standard
message-passing framework.

We propose to connect these two topics and build an explainable by-design
subgraph-enhanced GNN. We use a data-driven approach to learn small and
representative subgraphs that increase the expressive power for the downstream
task and that can be used as explanations.

##  2 Related Works

###  2.1 Explainability in GNNs

GNNs are generally seen as a black-box since they form decisions in high-
dimensional and nonlinear spaces. However, explainability techniques enable
humans to interpret the decision process of GNNs, discover potential sources
of error, find biases and limitations in the model and learn more about the
data and the task at hand. Most approaches in the literature are post-hoc
explainers differentiated according to the techniques used to explain a
trained GNN. In particular, there are gradient-based approaches [13],
perturbation-based approaches[5, 6, 14], decomposition methods [15], and
counterfactual explainers like [16] and GEM [17].

Of particular interest for this work is PGExplainer [6], which uses a small
network to parametrize the probability of each edge $\omega_{ij}$ of being
part of the explanatory subgraph, and sample from this distribution to obtain
the final explanation subgraph characterized by edges $e_{ij}$. The
optimization objective maximizes the mutual information between the masked
prediction and the prediction obtained from the original graph. In addition,
an element-wise entropy term encourages sparsity on $\omega$, and an $l1-$norm
forces small explanation subgraphs.

Fewer research efforts have been devoted to explainable-by-design GNNs.
ProtGNN [8] first learns some subgraph prototypes from the dataset. Then, they
compare the inputs to these prototypes in the latent space. Furthermore, they
introduce a subgraph sampling algorithm to highlight which component in the
input is most similar to each learned prototype. GIB [9], instead, uses a bi-
level optimization scheme to find the IB-subgraph, which is the most
informative yet compressed subgraph. This subgraph is the one that maximizes
the mutual information (MI) or shares the same properties with the original
input. GSAT [18], which builds upon GIB, injects stochasticity to the
attention weights to block the information from task-irrelevant graph
components while learning task-relevant subgraphs for interpretation.

###  2.2 Subgraph-enhanced graph neural networks

The study of the expressive power of GNNs has always been of central interest
to the community. Most GNNs operating via local neighbourhood aggregation are
as powerful as the Weisfeiler–Leman (1-WL) graph-isomorphism test [19]. It has
recently been shown that it is possible to create more expressive GNNs using
standard architectures that process several subgraphs of the input graph [11,
10]. In particular, ESAN [10] represents each graph as a bag of subgraphs
$\\{\mathcal{G}_{1},\dots,\mathcal{G}_{m}\\}$ chosen according to some
predefined policy, e.g., all graphs obtainable by removing one edge (edge
deleted strategy) or one node (node deleted strategy) from the original graph.
The encoder implements a module $L_{1}$ consisting of several message-passing
layer that process each subgraph independently, and then a second message-
passing module $L_{2}$ preprocesses the aggregation of the subgraphs working
as an information-sharing module for the subgraphs. For example, to compute
the embedding $H_{i}$ for the subgraph $\mathcal{G}_{i}$, they use the
following procedure:

$\text{H}_{i}=L_{1}\left(\mathcal{G}_{i}\right)+L_{2}\left(\sum_{j}\mathcal{G}_{j}\right)\,.$
(1)

Then, in the last layer of the encoder, a pooling operation aggregates the
node embeddings into subgraph embeddings with a global pooling operation.
Finally, a set learning module [20] aggregates the obtained subgraph
representations into a single one used in downstream tasks. The authors of
[10] name this general configuration DSS-GNN.

Recently k-OSANs [12] developed a data-driven policy for the subgraph
selection that improves the predictive performances compared to the simple
policies used in ESAN.

##  3 Proposed Framework

In this work, we consider an undirected and unweighted graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, where
$\mathcal{V}=\left\\{1,\ldots,n\right\\}$ is the set of node indexes, and
$\mathcal{E}{\color[rgb]{0.0,0.0,0.0}\subseteq}\left\\{(i,j)\;|\;i,j\in\mathcal{V}\right\\}$
is the set of edges connecting pairs of nodes. The entities and relationships
represented by nodes and edges depend on the data and the application. The
graph topology can be represented by the adjacency matrix
$\boldsymbol{\mathbf{A}}\in\left\\{0,1\right\\}^{n\times n}$. Other operators
matching the sparsity pattern of $\boldsymbol{\mathbf{A}}$, such as the graph
Laplacian, can be used to define the weighted connectivity of the graph.
Letting $\boldsymbol{\mathbf{D}}$ be one of these operators encoding the
topological information, a GNN layer is defined by (ignoring biases for
simplicity, but w.l.o.g.):

$\mathbf{H}=\phi\left(\mathbf{D}\mathbf{X}\mathbf{W}\right)\,,$ (2)

where $\mathbf{X}\in\mathbb{R}^{n\times d}$ is a matrix collecting all vertex
features row-wise, $\mathbf{W}\in\mathbb{R}^{d\times q}$ is a matrix of
trainable coefficients, and $\phi$ is an element-wise non-linear function
(e.g., for ReLU $\phi(s)=\max\left(0,s\right)$). A GNN can stack multiple
layers in the form of (2) to learn highly nonlinear node representations that
account for larger neighbourhoods on the graph. In this work, we are
interested in graph classification tasks. This task requires a global pooling
operation, such as mean, max, or min pooling, to build a graph representation
from the node representations produced by stacking the layers described in
Equation (2).

Figure 1: In the first step, we train the backbone of the SGNN, a GIN
classifier, using the original graphs. In the second step, we train the
explainer with the original graph and the backbone’s predictions. We then
create the new representation consisting of bags of explanation subgraphs.
Finally, in the third step, we train the whole SGNN framework, fine-tuning the
backbone using the explanatory subgraphs. Our model outputs the predicted
label and an explanation obtained by combining all the subgraphs used during
training.

Our goal is to develop a framework that jointly predicts the graph class and
the explanation masks, highlighting the parts of the graph that contribute the
most to the prediction. The key to our framework is the role played by the
subgraphs which are learned end-to-end, based on the optimization of the
classification loss. Firstly, they must improve the expressive power of the
base GNN classifier. Secondly, they must serve as explanation masks by
selecting the parts of the input that mostly contribute to determine the
correct class. We expect such subgraphs to be more informative than post-hoc
explanation masks since they are directly generated by the model to maximize
the classification performance.

We summarize our framework, which consists of three main steps highlighted in
Fig. 1. In the first step, we train the backbone of the SGNN, which is a
classifier that processes the original graph. In the second step we apply
PGExplainer with a minor modification to circumvent the “introduced evidence”
[21] issue due to the presence of soft masks. Specifically, we binarize the
soft weights to create “hard” explanation masks. Since “hard” explanations are
not differentiable, to allow the computation and the flowing of the gradients
in the backward pass, we apply a straight-through estimator [22]. These
modifications lead to design choices for the regularization terms that are
different from those originally used in PGExplainer. In particular, we apply
the regularization directly on the “hard” explanation mask $e$ instead of
$\omega$. Therefore we remove the element-wise entropy term and the $l1-$norm,
in this case, is equivalent to an $l0-$norm minimizing the number of edges
appearing in the explanation.

While PGExplainer generates a unique explanation subgraph, we are interested
in collecting several subgraphs to train the subgraph-based classifier.
Therefore, we devise two strategies to obtain a heterogeneous yet meaningful
set of subgraphs. The first is to use the reparametrization trick also in the
inference phase allowing the noise, which can be tuned as a hyperparameter, to
inject diversity that leads to multiple different subgraphs. In the second
approach we introduce a budget on the explanation. We adapt the binarization
threshold to obtain a set of explanations with $K$ edges, where $K$ ranges
from $5\%$ to $75\%$ of the initial number of edges in the original graph. An
example of these two strategies is presented in Figure 2, on the synthetic
dataset BA-2Motifs where the five nodes cycle is responsible for the
prediction.

Figure 2: The collection of subgraphs is represented as a heatmap. On the
left, there are the subgraphs obtained with the noise based strategy. On the
right, those obtained with the progressive top-K selection. The motif
responsible for the prediction is clearly visible in both cases, however, the
top-K approach selects less spurious edges.

In the third step, we insert the backbone GNN into the DSS-GNN model
introduced in ESAN in [10]. The backbone works as the encoder that processes
each explanation subgraph independently. A new message-passing module
preprocesses the aggregation of the explanations working as an information-
sharing module for the subgraphs. Finally, a new set learning module [20]
aggregates subgraph representations obtained after a global pooling operation
into a single one used in downstream tasks. Besides being used by the
information-sharing module, the subgraph obtained by aggregating all the
explanations is used to explain the model’s prediction. We would like to
stress the fact that this explanation, rather than being computed just from a
post-hoc explainer, is used to train the model.

##  4 Experimental Evaluation

We want to prove that our framework retains the same classification
performance of known SGNN architectures. However, our framework has the added
benefit of re-using the subgraphs involved in the computations as plausible
explanations for the model’s prediction. We used the same evaluation proposed
in ESAN [10] over the TUD repository datasets. This consist of computing the
average accuracy obtained by running a $10-$fold cross-validation. We selected
GIN [19] and DSS-GNN with the GIN backbone as baselines. We used the most
compact model consisting of 4 GIN layers with two linear layers of a hidden
size of 32. The batch size is also 32. We selected the top performing policies
from [10]: node deleted (ND) and edge deleted (ED). Furthermore, we consider
both their deterministic versions: one that uses the entire bag of subgraphs
(referred to as “1.0”) and the sampled versions, which keeps only $10\%$ of
the bag (referred to as “0.1”). Finally, we select our hyperparameters using
grid-search. We report the results in Table 1.

$\begin{array}[]{l|ccccc}\hline\cr\hline\cr\text{Model}&\text{MUTAG}&\text{NCI1}&\text{PROTEINS}&\text{IMDB-
BINARY}&\text{IMDB-MULTI}\\\ \hline\cr\text{GIN}&89.4\pm 5.6&82.7\pm
1.7&76.2\pm 2.8&75.1\pm 5.1&52.3\pm 2.8\\\ \text{DSS ND }1.0&89.0\pm
4.4&83.0\pm 2.0&76.5\pm 3.7&75.4\pm 3.5&52.3\pm 2.9\\\ \text{DSS ND
}0.1&89.9\pm 2.8&\textbf{83.1}\pm 2.0&\textbf{76.6}\pm 4.7&75.7\pm 2.3&52.5\pm
2.2\\\ \text{DSS ED }1.0&89.9\pm 4.3&82.7\pm 1.7&74.0\pm 4.4&\textbf{76.7}\pm
3.6&52.6\pm 2.6\\\ \text{DSS ED }0.1&89.9\pm 5.5&82.1\pm 2.0&75.1\pm
4.7&76.4\pm 2.7&52.6\pm 2.8\\\ \text{OURS}&\textbf{92.1}\pm 4.3&82.9\pm
2.3&76.3\pm 4.4&76.0\pm 3.1&\textbf{53.5}\pm 3.8\\\
\hline\cr\hline\cr\end{array}$ Table 1: TUDataset classification accuracy. We
report the mean and the standard deviation obtained by running a $10-$fold
cross-validation.

We can observe that subgraph-enhanced GNNs achieve slightly better accuracies.
We did not notice any empirical difference when using a subset of the bag of
subgraphs instead of the complete one. Furthermore, as noted in [10], despite
the model being no longer invariant, the sample bag of subgraphs still
increases the expressive power of the base encoder. Finally, the results
obtained using explanatory subgraphs instead of randomly generated ones are
comparable. However, we have the advantage that our bag of subgraphs forms an
explanation helping human experts to extract new knowledge from the model.

##  5 Conclusion

We introduced a framework to bridge subgraph-enhanced GNNs with explainability
techniques. First, we pre-train a base GNN such that we can train a post-hoc
explainer modified to satisfy the need for multiple and discrete explanation
masks. Next, the explainer generates a bag of explanation subgraphs. Then, we
resume the training of the base model, including equivariant layers, to
process the new representation in the SGNN fashion. Finally, we combine the
subgraphs in a unique heatmap highlighting the most relevant parts responsible
for the prediction and the prediction itself. This additional information
allows human experts to interpret and extract knowledge from the model. The
experimental evaluation showed that our model achieves comparable
classification performance to state-of-the-art models while providing, at the
same time, an explanation.

## References

  * [1] Davide Bacciu, Federico Errica, Alessio Micheli, and Marco Podda.  A gentle introduction to deep learning for graphs.  Neural Networks, 129:203–221, 2020.
  * [2] Marco Gori, Gabriele Monfardini, and Franco Scarselli.  A new model for learning in graph domains.  In Proc. 2005 IEEE International Joint Conference on Neural Networks (IJCNN), volume 2, pages 729–734. IEEE, 2005.
  * [3] Alessio Micheli.  Neural network for graphs: A contextual constructive approach.  IEEE Transactions on Neural Networks, 20(3):498–511, 2009.
  * [4] Max W. Shen.  Trust in ai: Interpretability is not necessary or sufficient, while black-box interaction is necessary and sufficient.  ArXiv, abs/2202.05302, 2022.
  * [5] Rex Ying, Dylan Bourgeois, Jiaxuan You, Marinka Zitnik, and Jure Leskovec.  Gnnexplainer: Generating explanations for graph neural networks.  Advances in Neural Information Processing Systems, 32:9240, 2019\.
  * [6] Dongsheng Luo, Wei Cheng, Dongkuan Xu, Wenchao Yu, Bo Zong, Haifeng Chen, and Xiang Zhang.  Parameterized explainer for graph neural network.  Advances in Neural Information Processing Systems, 33, 2020.
  * [7] Indro Spinelli, Simone Scardapane, and Aurelio Uncini.  A meta-learning approach for training explainable graph neural networks.  IEEE Transactions on Neural Networks and Learning Systems, pages 1–9, 2022.
  * [8] Zaixi Zhang, Qi Liu, Hao Wang, Chengqiang Lu, and Cheekong Lee.  Protgnn: Towards self-explaining graph neural networks.  CoRR, abs/2112.00911, 2021.
  * [9] Junchi Yu, Tingyang Xu, Yu Rong, Yatao Bian, Junzhou Huang, and Ran He.  Graph information bottleneck for subgraph recognition.  In ICLR, 2021.
  * [10] Beatrice Bevilacqua, Fabrizio Frasca, Derek Lim, Balasubramaniam Srinivasan, Chen Cai, Gopinath Balamurugan, Michael M. Bronstein, and Haggai Maron.  Equivariant subgraph aggregation networks.  CoRR, abs/2110.02910, 2021.
  * [11] Leonardo Cotta, Christopher Morris, and Bruno Ribeiro.  Reconstruction for powerful graph representations.  CoRR, abs/2110.00577, 2021.
  * [12] Chendi Qian, Gaurav Rattan, Floris Geerts, Christopher Morris, and Mathias Niepert.  Ordered subgraph aggregation networks, 2022.
  * [13] Phillip E. Pope, Soheil Kolouri, Mohammad Rostami, Charles E. Martin, and Heiko Hoffmann.  Explainability methods for graph convolutional neural networks.  In 2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 10764–10773, 2019.
  * [14] Michael Sejr Schlichtkrull, Nicola De Cao, and Ivan Titov.  Interpreting graph neural networks for NLP with differentiable edge masking.  CoRR, abs/2010.00577, 2020.
  * [15] Thomas Schnake, Oliver Eberle, Jonas Lederer, Shinichi Nakajima, Kristof T Schutt, Klaus-Robert Mueller, and Gregoire Montavon.  Higher-order explanations of graph neural networks via relevant walks.  IEEE Transactions on Pattern Analysis and Machine Intelligence, pages 1–1, 2021.
  * [16] Ana Lucic, Maartje ter Hoeve, Gabriele Tolomei, Maarten de Rijke, and Fabrizio Silvestri.  Cf-gnnexplainer: Counterfactual explanations for graph neural networks.  arXiv preprint arXiv:2102.03322, 2021.
  * [17] Wanyu Lin, Hao Lan, and Baochun Li.  Generative causal explanations for graph neural networks.  CoRR, abs/2104.06643, 2021.
  * [18] Siqi Miao, Miaoyuan Liu, and Pan Li.  Interpretable and generalizable graph learning via stochastic attention mechanism.  CoRR, abs/2201.12987, 2022.
  * [19] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  In International Conference on Learning Representations, 2019.
  * [20] Haggai Maron, Or Litany, Gal Chechik, and Ethan Fetaya.  On learning sets of symmetric elements.  In Proceedings of the 37th International Conference on Machine Learning, ICML’20. JMLR.org, 2020.
  * [21] Piotr Dabkowski and Yarin Gal.  Real time image saliency for black box classifiers.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [22] Geoffrey Hinton, Li Deng, Dong Yu, George E. Dahl, Abdel-rahman Mohamed, Navdeep Jaitly, Andrew Senior, Vincent Vanhoucke, Patrick Nguyen, Tara N. Sainath, and Brian Kingsbury.  Deep neural networks for acoustic modeling in speech recognition: The shared views of four research groups.  IEEE Signal Processing Magazine, 29(6):82–97, 2012.


# Microtubules as electron-based topological insulators

Varsha Subramanyan<EMAIL_ADDRESS>Dept. of Physics, University of
Illinois, Urbana-Champaign,USA  Kay L. Kirkpatrick  Dept. of Mathematics,
University of Illinois, Urbana-Champaign,USA  Dept. of Physics, University of
Illinois, Urbana-Champaign,USA  Saraswathi Vishveshwara  Molecular Biophysics
Unit, Indian Institute of Science, Bengaluru, India  Smitha Vishveshwara
<EMAIL_ADDRESS>Dept. of Physics, University of Illinois, Urbana-
Champaign,USA

###### Abstract

The microtubule is a cylindrical biological polymer that plays key roles in
cellular structure, transport, and signalling. In this work, based on studies
of electronic properties of polyacetelene and mechanical properties of
microtubules themselvesProdan and Prodan (2009), we explore the possibility
that microtubules could act as topological insulators that are gapped to
electronic excitations in the bulk but possess robust electronic bounds states
at the tube ends. Through analyses of structural and electronic properties, we
model the microtubule as a cylindrical stack of Su-Schrieffer-Heeger chains
(originally proposed in the context of polyacetylene) describing electron
hopping between the underlying dimerized tubulin lattice sites. We postulate
that the microtubule is mostly uniform, dominated purely by GDP-bound dimers,
and is capped by a disordered regime due to the presence of GTP-bound dimers
as well. In the uniform region, we identify the electron hopping parameter
regime in which the microtubule is a topological insulator. We then show the
manner in which these topological features remain robust when the hopping
parameters are disordered. We briefly mention possible biological implications
for these microtubules to possess topologically robust electronic bound
states.

Organic polymers and biomolecules have served as inspiration for discovering
new condensed matter phenomena, with path-breaking insights developed for both
biological systems and unique physical modelsSpakowitz (2019); ALVAREZ-ESTRADA
and CALVO (2002). The Su-Schrieffer-Heeger (SSH) model proposed in the 1980’s
for polyacetylene serves even today as the simple paradigm for topological
phases of matter and associated charge fractionalizationSu  _et al._ (1980).
Such models have then led to opening entire fields of study, bringing
theoretical and experimental advances in our understanding of materialsAsbóth
_et al._ (2016). On the biological side, the fundamental units of the model,
involving dimerization and sublattice symmetry, are the building blocks for a
range of complex macromolecules. Here, we investigate microtubules as such a
highly promising instance for extending topological behavior in rich and
diverse ways.

The microtubule is a crucial protein complex that forms an integral part of
the cellular cytoskeleton and consists of dimerized units arranged in a
helical structureWade (2009); Amos (2011). In Ref. Prodan and Prodan, 2009,
the authors target the mechanical properties of these dimer units embedded in
a manifold in the uniform (non-disordered) case and model the underlying
interactions as harmonic coupling. The resulting dynamical equations show
similarities with the Hamiltonian of a topological insulator. Strikingly, the
similarity translates to the mechanical system being topological in that its
phases can be characterized by certain invariants and topologically non-
trivial regimes in parameter space possess a class of solutions reflecting
signature end modes rendered robust by a spectral gap. As with
metamaterialsXin  _et al._ (2020), the authors associate these modes with
robust boundary phonon modes that play an important role in the growth and
decay phases of microtubules. We do the same in a corresponding electronic
system, motivated by the healthy interplay between metamaterial and electronic
physics. The crucial difference is that here the systems are one and the
same–we expect that the predicted phonon behavior in microtubules would imply
corresponding analogues in its electronic properties.

Figure 1: The cylindrical microtubule is shown here in different phases and
configurations. It is composed mainly of GDP-bound dimers (shades of green).
Figure a) shows a microtubule of zero pitch, while Figure b) shows a
microtubule having pitch equal to one dimer’s length. Figures c) and d) show
the microtubule in two possible phases. While a) and b) are uniform, in c),
the “cap” of the microtubule has dimers of both the GDP and GTP (blue) kind,
leading to increase in its length when more dimers from the surrounding medium
attach onto the cap. This non-uniformity prevents the peeling off or reduction
in microtubule length, depicted in Figure d). Figure e) shows the lattice
model of the microtubule with associated hopping strengths indicated. In
Figure f) the $\alpha$ and $\beta$ tubulin monomers are indicated in each type
of dimer.

The attribution of topological properties to electronic systems has important
consequences Bernevig (2013). For instance, fractional charge polarization,
originally analyzed in polyacetelene Meier  _et al._ (2016)and a zero bias
peak in the differential conductivity are defining signatures of topological
phases. These quantities are also robust to small amounts of disorder
permeating the system due to the presence of a finite energy gap. However, it
is also known that stronger disorder strengths mediate phase transitions
between topological and trivial phases.

In this work, we show that within a regime of plausible parameter windows,
microtubules can possess a topological phase determined by their electronic
band structure, alongside the mechanical phase posited in Ref. Prodan and
Prodan, 2009. On consideration of electronic processes inherent to the
microtubule, we model it as a cylindrical stack of SSH chains. We explore the
properties of such quasi one-dimensional weak topological insulators and their
response to disorder. We then discuss how these electronic phases map back to
biological properties, possibly showing radically different kinds of
conduction properties of microtubules stemming from the presence of
topologically robust bound states.

Biophysics, structure, and electronic properties:  Microtubules play important
roles in maintaining structural integrity of the cell, material transport
within it, in initiating cellular mitosis, and in facilitating intracellular
signallingAvila (1992). They are prominent structures in eukaryotic cells,
especially crucial to the morphological make-up of neurons. The microtubule
has long been a popular candidate for the exploration of quantum mechanical
effects manifesting in biological systemsHameroff and Watt (1982); Hameroff
and Penrose (1996). It has often been argued that finite temperature and the
environment associated with living cells would lead to decoherence and
suppression of such effects. However, several recent studies provide evidence
that such effects may not be ignored in several biological systems including
microtubulesMarais  _et al._ (2018); Craddock  _et al._ (2014); Mavromatos and
Nanopoulos (1998). Here, we base our fundamental assumptions on such evidence.

Structurally, the microtubule is a helical cylinder composed of tubulin
molecules as the fundamental unit. These units are dimers consisting the
homologous monomers $\alpha-$tubulin and $\beta-$tubulin. Typically, the
cylinders have a 25 nm outer diameter and 15 nm inner diameter. The dimers
themselves have dimensions of 8 nm $\times$ 5 nm $\times$ 5 nmTuszynski  _et
al._ (2018). The dimers are attached helically around the cylinder giving rise
to protofilaments, usually 13 in number around the circumference. The
cylindrical lattice thus formed is highly ordered and is most commonly seen in
two different configurations - with and without a prominent seam,
corresponding to differences in helical pitch, as seen in Fig 1. As has been
done with regards to mechanical properties in Ref. Prodan and Prodan, 2009,
these cylinders can thus be modeled as extensions of the SSH system where the
chain is wrapped around in a cylindrical configuration.

The microtubule shows two distinct phases when it undergoes “dynamic
instability” - a growth or polymerization phase, where dimers from the
surrounding solution attach themselves to the filaments of the microtubule,
and a “peeling off” or de-polymerization or catastrophe phase, where the
filaments of the microtubule curve outward longitudinally and lose dimers to
the solution, decreasing the effective length of the tubeBrouhard (2015);
Horio and Murata (2014). The microtubule has been observed to stochastically
switch between the two phases. The dimers themselves can be of two possible
types distinguished by the molecule attached at an exchange site - GDP or GTP.
During the growth phase, dimers of both kinds attach themselves to the
microtubule to form longer filaments. This non-uniform segment at the end of
tube is called the cap. As the tube grows, the GTP hydrolyzes to GDP after
attachment to form uniform GDP attached extended regions. The GDP-attached
filaments are longitudinally curved, and it has been proposed that this
curvature that results in the “peeling” off when the energetics overcome the
stability provided by the capBurbank and Mitchison (2006); Brun  _et al._
(2009). Based on this structure of the microtubule, in our model, we
distinguish two distinct regimes—one in the uniform GDP-based region, which we
model as consisting of no disorder, and the other in the cap region, which we
model as having a disorder distribution in the inter-tubulin charge hopping
amplitudes.

Turning to the electronic properties, ample computational and observational
studies indicate that the dimerized structure of the microtubule affect not
only vibrational degrees of freedom, but also the intersite electronic hopping
within the underlying latticeEakins  _et al._ (2021a); Tuszynski (2019);
Minoura and Muto (2006). The conduction electron concentration on a
microtubule has been estimated to be $n\sim 10^{19}cm^{-3}$, which is higher
than that of semiconductors, but significantly lower than that of metals. Even
lower estimates of its conductivity from these numbers (at physiological
temperatures) suggest a conductivityTuszynski  _et al._ (2018) that lies
between $0.04\Omega^{-1}m^{-1}$ and $10^{5}\Omega^{-1}m^{-1}$. The microtubule
has a large net negative charge per tubulin dimer that is mostly concentrated
on the ”C-terminus” of each dimer. The uncompensated negative charge is
balanced in solution by counterions that screen the charge, giving rise to a
dielectric polarization around the surface of the microtubule. This
observation has led to several models of conductivity in terms of tight-
binding of electrons to the lattice site and effective charge hopping from
site to site, like in the case of organic semiconductorsRühle  _et al._
(2011); Kirkpatrick  _et al._ (2013). Additionally, since the surrounding
medium of the microtubule is rich in positively charged ionsTuszynski (2019);
Eakins  _et al._ (2021b), the highly negatively charged microtubule might also
exhibit conduction by proton hopping mechanisms, leading to higher
conductivity, a phenomenon that plays a secondary role in this work. These
descriptions for conductivity form the starting point of our analysis; we
model electronic properties in terms of electrons hopping between nearest
neighbor lattice sites provided by the underlying $\alpha$\- and
$\beta$-tubulin units.

Effective Quasi-1D model: Having established this tight-binding description
for electron hopping between the $\alpha$ and $\beta$ sites of the
microtubule, we first analyze the uniform region. The resultant cylindrical
lattice and hopping strengths are as depicted in Fig. 1. The system consists
of SSH chains along the length of the tube having alternating bond strengths
$t_{1}$ and $t_{2}$. As shown in the figure, along the circumferential
direction, the units are connected via bonds $t_{\perp,1}$ and $t_{\perp,2}$.
Depending on the pitch of the helix, there will be conformational changes in
the helix, which can be reflected in the relative strength between
$t_{\perp,1}$ and $t_{\perp,2}$. The pitch will also be reflected in the
allowed quantized values of $k_{y}$. Here, due to the helical nature of the
tubule, we assume that the staggering of the chains results in these inter-
chain hoppings dominating. Thus, the predominant hopping terms respect sub-
lattice symmetry in that hoppings are only between sub-lattices $A$ consisting
of $\alpha$ sites and $B$ consisting of $\beta$ sites, but not within each
sublattice. The resultant Hamiltonian thus takes the generic form:

$\displaystyle H$
$\displaystyle=\sum_{ij}\begin{pmatrix}c_{A,ij}^{\dagger}&c_{B,ij}^{\dagger}\end{pmatrix}\mathbb{H}\begin{pmatrix}c_{A,ij}\\\
c_{B,ij}\end{pmatrix}$ (1)

where $\mathbb{H}$ obeys chiral symmetry such that $S\mathbb{H}S=-\mathbb{H}$
with $S=\sum_{ij}(c^{\dagger}_{A,ij}c_{A,ij}-c^{\dagger}_{B,ij}c_{B,ij})$.
This Hamiltonian, emerging from the dominant hoppings, respects the symmetries
of quasi one dimensional topological systems belonging to class AIII,
classified as per the usual Altland-Zirnbauer formalism that characterizes
topological systems according to their symmetriesSchnyder  _et al._ (2009).

A qualitative picture of how an ’end-mode’ can form along the circumference at
the two ends of a segment of a uniform region is as follows. In the absence of
interchain hopping, the paradigm SSH chains running along the length of the
tube naturally support end bound states. These bound states correspond to
robust mid-gap states in a spectrum symmetric about the Fermi energy having a
gap of magnitude $|t_{1}-t_{2}|$. Upon inter-chain coupling, these states too
couple to form the end-mode along the rimDelplace  _et al._ (2011); Zhu  _et
al._ (2019); the robustness of similar states has been discussed in the
context of coupled ladder systems Li  _et al._ (2017); Hügel and Paredes
(2014); Padavić  _et al._ (2018). Figure.2 b) and c) plot such typical end-
mode eigenstates in a regime that we have determined to be topological. For
strong enough interchain coupling, the gap in the spectrum closes and the
system undergoes a transition into a phase where no such robust end-mode
exists.

To analyze these features more rigorously, we express the Hamiltonian for a
segment of this lattice in momentum space. To first order, we assume that the
segment is long enough to take the thermodynamic limit. We assume a unit
lattice constant for each bond. The resultant Hamiltonian takes the form

$\displaystyle H$
$\displaystyle=\begin{pmatrix}c_{A,\vec{k}}^{\dagger}&c_{B,\vec{k}}^{\dagger}\end{pmatrix}\begin{pmatrix}0&q(\vec{k})\\\
q^{\dagger}(\vec{k})&0\end{pmatrix}\begin{pmatrix}c_{A,\vec{k}}\\\
c_{B,\vec{k}}\end{pmatrix}$ (2)

where $q(\bar{k})=u+ve^{-ik_{x}}$, $\bar{k}=(k_{x},k_{y})$,
$u=t_{1}+2t_{\perp,2}\cos{k_{y}}$ and $v=2t_{\perp,1}\cos{ky}+t_{2}$. The
energy band structure is given by

$\displaystyle E_{\pm}$ $\displaystyle=\pm\sqrt{u^{2}+v^{2}+2uv\cos{k_{x}}}.$
(3)

In order to chart out the regimes where the system is topological in that a
finite-sized tube would have end modes present, we can use the standard
topological analysis of the band structure. In the thermodynamic limit, or
under periodic boundary conditions in both directions, the band structure is
gapped for a range of parameter regimes, as shown in the phase diagram in Fig.
2. The gapped phases may further be classified as topological or trivial, as
depicted in the phase diagram. Crucially, the topological phase is
characterized by a topological invariant or winding number defined as

$\displaystyle\nu_{x}=-\frac{i}{2\pi}\int dk_{x}Tr[q^{-1}\partial_{k_{x}}q].$
(4)

In this system, $\nu$ takes the value 1 for $|u|<|v|$ (topological) or 0 for
$|u|>|v|$ (trivial)Asbóth  _et al._ (2016). The system is gapless when
$|u|=|v|$.

Figure 2: Energy bands and phase diagram: Figure a) shows the phase diagram
for the lattice model under the thermodynamic limit, as a function of hopping
strengths $t_{2}$, $t_{\perp,1}$ and $t_{\perp,2}$ all scaled to be in units
of $t_{1}$.The gapless and gapped phases are indicated, with the gapped phases
being identified by their associated topological invariant. Figures b) and c)
show one of the possible zero energy edge states for a cylindrical strip of
$N_{y}=13$ chains with 50 dimers along the length of the microtubule in the
topological phase. Figure b) is for a helical cylinder of pitch zero, while c)
is for a helical cylinder with unit pitch. Figures d) and e) demonstrate the
changes in the band structure for strip of pitch zero (Figure d)) and pitch
one (Figure e)) when the parameter $t_{2}$ is tuned from 0 to 3. Here,
$t_{\perp,1}=0.2$ and $t_{\perp,2}=0.1$. The figure shows an initial
configuration with a trivial gap become gapless for a range of values before
again opening up a topological gap. The zero energy modes are also seen in the
topological gap.

The band structure of a finite cylindrical strip is shown in Fig. 2 d) and e),
for two different boundary conditions of the microtubule. In the figures, as
parameter $t_{2}$ is increased, it is seen that that system initially has a
trivial gap, and then passes through a regime of gaplessness before opening up
a topological gap. The topological phase is distinguished by the presence of
zero energy (compared to the Fermi energy) edge states. For a cylindrical
strip of $N_{y}$ chains, there are $2N_{y}$ zero energy modes in the gap, in
the absence of sub-lattice symmetry breaking terms. In our quasi 1D model,
this situation corresponds to a topological invariant $\nu_{x}=1$ (higher
integers are usually associated with long-range coupling in such conjugated
systems). These edge states can then couple with other systems, hybridize with
each other for short enough tube lengths, or interact with environmental
modes. These effects give rise to a host of phenomena that have been
extensively studied in the context of Hermitian as well as non-Hermitian SSH
chains and would be of direct relevance to microtubulesLieu  _et al._ (2020);
McGinley and Cooper (2019); Campos Venuti  _et al._ (2017). Moreover, while
the staggered helical nature of the microtubule allows the reasonable
assumption of low hoppings of the kind A-A or B-B or non-uniform on-site
chemical potential, their presence explicitly breaks the aforementioned
sublattice symmetry of the system. Its effects may be understood in terms of
the effects of extended nest-nearest-neighbour coupling in SSH chainsPérez-
González  _et al._ (2019). For small values of such symmetry breaking terms
(with respect to the gap), the spectrum is shifted or stretched or both, thus
lifting some or all the edge states from zero energy. However, they retain
their topological significance despite due to the presence of the gap. Very
high values, on the other hand, can close the gap and cause these edge states
to become absent.

Disorder and Dynamic Instability: While the above analysis is most relevant to
the uniform regime corresponding to the peeling-off segment, by the cap
theory, the growth phase of the microtubule sets in when the ends are not
uniformly composed of GDP dimers, but also contain interspersed GTP dimers.
Further, in some diseased states, there could also be isoforms of tubulin
present, which would further contribute to the non-uniformity of the capVemu
_et al._ (2017). This non-uniformity can be reasonably expected to affect the
electron hopping strength. Moreover, given that the GTP molecules are in
metastable stages of hydrolysis enroute to conversion to GDP, we expect a
continuous range of non-uniformity, as opposed to a binary distribution. We
model the non-uniformity in the cap as disorder in the hopping parameters, and
characterize its effects on topological properties via changes made to the
winding number.

Since the cap does not cover the entire length of the microtubule, we assume
disorder only in part of the tube. As a specific proof-of-principle instance,
we model the presence of disorder in half the tube. In order to determine
whether or not the tube is in the topological phase, we can no longer employ
the momentum basis form for deriving the invariant of Eq.4 since disorder
breaks translational invariance. We may however employ the real space version
typically used for disordered systems Prodan (2011); Claes and Hughes (2020);
Mondragon-Shem  _et al._ (2014). To briefly provide the recipe, in the basis
of eigenstates of $S$, the Hamiltonian can be written in the off-diagonal form
of Eq.2. One then defines a ’flat-band’ version $Q$ of this real-space
Hamiltonian as $Q=P_{+}-P_{-}$, where $P_{\pm}$ are projectors onto
eigenspaces with positive and negative eigenvalues. The invariant is given by

$\displaystyle\nu=-\frac{1}{N_{x}}Tr(\tilde{Q}^{\dagger}[X,\tilde{Q}]),\textnormal{
}Q=\begin{pmatrix}0&\tilde{Q}\\\ \tilde{Q}^{\dagger}&0\end{pmatrix}.$ (5)

Here $N_{x}$ is the number of unit cells along the lattice and $X$ is the
position operator in the lattice.

Figure 3: Topological invariant $\nu_{x}$ plotting against disorder strength W
in $t_{1}$ and strength W/2 in $t_{2}$ for $N_{y}$=8 in two possible pitch
configurations. The blue points correspond to the real space topological
invariant for a microtubule of pitch zero, and the red curve is a rough fit of
the data points. The initial quantization is briefly retained until $W\sim
3.5$. The system briefly stabilizes at $\nu_{x}=3/8$ before the gap closes for
stronger disorder.

Our model of the microtubule corresponds to the AIII class of topological
insulators. We effectively evaluate the topological invariant per unit width
$\nu_{x}=\nu/N_{y}$ where $N_{y}$ is the number of protofilaments. Thus,
$\nu_{x}$ is quantized in units of $1/N_{y}$. This quantized behaviour has
been established both analytically and numerically, even in the presence of
strong disorderClaes and Hughes (2020). In our system, we apply disorder to
the hopping strengths $t_{1}$ and $t_{2}$ individually. For each of these
hopping strengths (indicated generally as $h_{i}$), disorder is applied
through the equation $h_{i}=h_{i0}+W_{i}u$ where $u$ is a uniform random
variable sampled between $[-0.5,0.5]$ and $W_{i}$ is the strength of the
disorder. The initial parameters $h_{i0}$ are picked so that the microtubule
starts at a topological phase. We also pick $W_{1}=W$ for $t_{1}$ and
$W_{2}=W/2$ for $t_{2}$. These parameters represent a typical sample which
demonstrates the plausible evolution of topological features.

The results thus obtained are shown in Fig.3. A few caveats are in order.The
trend clearly shows a nearly quantized value of unity for the invariant at low
disorder that drops towards zero upon increasing disorder. The curve can
stabilize at other values of $\nu_{x}$ that are multiples of $1/N_{y}$.
Fluctuations around the allowed values are a result of the finite size of the
system chosen and they diminish for longer systemsClaes and Hughes (2020).
Given the computational power involved, our strips were confined to $N_{y}=8$.
We expect the results to hold for the more realistic number $N_{y}\sim 13$.
While the form of the invariant assumes sublattice symmetry, we expect the
topological feature of a robust edge mode to remain for small deviations from
this symmetry condition. Most importantly, the simulation is a demonstration
of the salient properties in the presence of disorder. Our simulations show
that the quantization is retained for values of disorder as high as $W=3.5$.
While this treatment considers the segment as a whole, with disorder in half
the region to represent the physical situation, a more detailed analysis would
show the exact physical location of the end modes as possibly shifting into
the interior. Most prominently, we may reasonably expect that for regions of
lower disorder, the presence of the edge mode would still be prevalent.

Discussion and Outlook: Here, we have given a proof-of-concept presentation of
how topological properties of tight-binding lattices can provide biological
microtubules with unique electronic bound state structure hitherto unexplored.
In particular, we associate the decay phase of the microtubule with the
presence of an electronic topological edge mode, while we associate the growth
phase with the disordered case that could possibly lead to the vanishing of
the edge mode. Several variations of course would come into play with regards
to the rudimentary model presented here. To mention a few, in this work, we
assumed sublattice symmetry, which ultimately gives rise to the mid-gap states
centered at the Fermi energy. We expect that deviations either due to on-site
potential variations or other couplings within each of the $A$ and $B$
sublattices would still preserve the edge modes for small enough values. On
the other hand, for large enough deviations, these edge states would
disappear; perhaps the asymmetry of the physical properties in the two ends of
the tubes could emerge from related differences in the chemical composition.
Additionally, one could adapt these calculations to extremely short
microtubules as they start to grow out of the centrosome, and perhaps to
extremely long microtubules that push an embryonic neuronal cell membrane out
to grow towards future synapses. Two other important factors to consider are
coupling to the environment and the non-equilibrium nature of tube growth and
dynamic instability.

With regards to disorder, while our model is phenomenological in nature,
further studies would connect the disorder distribution to physical
conditions, such as hydrolysis states from GTP to GDP, temperature effects,
and presence of other biochemicals. Moreover, here, while we studied uniform
versus disordered behavior for entire segments, it is possible to have more
robust bound states at the interface of such segments. It is postulated that
when isotypes of tubulin are present in diseased microtubules, their
functioning is only marginally altered, with such defects only becoming
apparent in processes that are highly sensitive to microtubule
dynamicsGadadhar  _et al._ (2017). This view additionally supports the
presence of topologically robust modes in the microtubule. Another important
consideration is the previous study on mechanical degrees of freedom
demonstrating that the uniform region is unstable due to a topological phonon
mode, in contrast to the electronic modes studied here. We expect both modes
to follow similar trends in the presence of disorder. A further intriguing
possibility would entail electron-phonon coupling.

Finally, the existence of robust edge bound state modes could have significant
biological implications. In principle, the bound state could act as an
acceptor or donor of charge, perhaps differentially influencing the more
negatively charged GTP compared with the GDP molecule for the processes
discussed here. The bound state could also be affected by microtubule
associated proteins, such as as kinesin and dynein. Microtubules are key to
intracellular long-range transport, for instance, quite strikingly in neural
transmission. In addition, if edge modes of microtubules provide a localized
state for electronic occupation, then they may extend into various parts of
the axon terminal, and dynamically change the open-or-closed states of
voltage-gated ion channels, contributing to pre-synaptic membrane activity.
Conversely, proteins that have an overall positive charge or neurite regions
that have a concentration of positive charge (such as an action potential’s
influx of cations) might attract the C-terminus of a microtubule if it is in
the topological regime for hosting negative charge at its end. Regardless of
the extent to which microtubules are involved with intra-neuronal activity
besides transport Hameroff and Watt (1982); Hameroff and Penrose (1996), it is
likely that sub-neuronal computation is an important part of brain function,
and microtubules may play a significant role in that. As in condensed matter
systems, conduction properties of microtubules hinge greatly on their
structure, chemical composition and environment; their possible behavior as
topological insulators may bear radical insights into the quantum nature of
biophysical systems.

Acknowledgements: We thank Camelia Prodan and Emil Prodan for sharing
beautiful insights both in their work and through conversations. We also thank
Kesav Krishnan for useful conversations. We acknowledge the support of the
National Science Foundation under grants DMR-2004825 (V. S. and Sm.V.) and
CAREER DMS-1254791 (K. L. K.). Sa.V. is an Honorary Scientist of the National
Academy of Sciences, India(NASI).

## References

  * Prodan and Prodan (2009) E. Prodan and C. Prodan, Phys. Rev. Lett. 103, 248101 (2009).
  * Spakowitz (2019) A. J. Spakowitz, The Journal of Chemical Physics 151, 230902 (2019), https://doi.org/10.1063/1.5126852 .
  * ALVAREZ-ESTRADA and CALVO (2002) R. F. ALVAREZ-ESTRADA and G. F. CALVO, Molecular Physics 100, 2957 (2002), https://doi.org/10.1080/00268970210121623 .
  * Su  _et al._ (1980) W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. B 22, 2099 (1980).
  * Asbóth  _et al._ (2016) J. K. Asbóth, L. Oroszlány, and A. Pályi, Lecture Notes in Physics (2016), 10.1007/978-3-319-25607-8.
  * Wade (2009) R. H. Wade, Molecular Biotechnology 43, 177 (2009).
  * Amos (2011) L. A. Amos, Semin Cell Dev Biol 22, 916 (2011).
  * Xin  _et al._ (2020) L. Xin, Y. Siyuan, L. Harry, L. Minghui, and C. Yanfeng, Current Opinion in Solid State and Materials Science 24, 100853 (2020).
  * Bernevig (2013) B. A. Bernevig,  _Topological Insulators and Topological Superconductors_ (Princeton University Press, 2013).
  * Meier  _et al._ (2016) E. J. Meier, F. A. An, and B. Gadway, Nature Communications 7, 13986 (2016).
  * Avila (1992) J. Avila, Life Sci 50, 327 (1992).
  * Hameroff and Watt (1982) S. R. Hameroff and R. C. Watt, Journal of Theoretical Biology 98, 549 (1982).
  * Hameroff and Penrose (1996) S. Hameroff and R. Penrose, Mathematics and Computers in Simulation 40, 453 (1996).
  * Marais  _et al._ (2018) A. Marais, B. Adams, A. K. Ringsmuth, M. Ferretti, J. M. Gruber, R. Hendrikx, M. Schuld, S. L. Smith, I. Sinayskiy, T. P. J. Krüger, F. Petruccione, and R. van Grondelle, Journal of The Royal Society Interface 15, 20180640 (2018).
  * Craddock  _et al._ (2014) T. J. A. Craddock, D. Friesen, J. Mane, S. Hameroff, and J. A. Tuszynski, Journal of The Royal Society Interface 11, 20140677 (2014).
  * Mavromatos and Nanopoulos (1998) N. E. Mavromatos and D. V. Nanopoulos, International Journal of Modern Physics B 12, 517 (1998), https://doi.org/10.1142/S0217979298000326 .
  * Tuszynski  _et al._ (2018) J. A. Tuszynski, A. Priel, J. A. Brown, H. F. Cantiello, and J. M. Dixon, Nano and Molecular Electronics Handbook (2018).
  * Brouhard (2015) G. J. Brouhard, Molecular biology of the cell 26, 1207 (2015), 25823928[pmid].
  * Horio and Murata (2014) T. Horio and T. Murata, Frontiers in Plant Science 5, 511 (2014).
  * Burbank and Mitchison (2006) K. S. Burbank and T. J. Mitchison, Curr Biol 16, R516 (2006).
  * Brun  _et al._ (2009) L. Brun, B. Rupp, J. J. Ward, and F. Nédélec, Proceedings of the National Academy of Sciences 106, 21173 (2009), https://www.pnas.org/content/106/50/21173.full.pdf .
  * Eakins  _et al._ (2021a) B. B. Eakins, S. D. Patel, A. P. Kalra, V. Rezania, K. Shankar, and J. A. Tuszynski, Frontiers in Molecular Biosciences 8, 150 (2021a).
  * Tuszynski (2019) J. A. Tuszynski, “The bioelectric circuitry of the cell,” in  _Brain and Human Body Modeling: Computational Human Modeling at EMBC 2018_, edited by S. Makarov, M. Horner, and G. Noetscher (Springer International Publishing, Cham, 2019) pp. 195–208.
  * Minoura and Muto (2006) I. Minoura and E. Muto, Biophysical journal 90, 3739 (2006), 16500962[pmid].
  * Rühle  _et al._ (2011) V. Rühle, A. Lukyanov, F. May, M. Schrader, T. Vehoff, J. Kirkpatrick, B. Baumeier, and D. Andrienko, Journal of Chemical Theory and Computation 7, 3335 (2011).
  * Kirkpatrick  _et al._ (2013) K. Kirkpatrick, E. Lenzmann, and G. Staffilani, Communications in Mathematical Physics 317, 563 (2013).
  * Eakins  _et al._ (2021b) B. B. Eakins, S. D. Patel, A. P. Kalra, V. Rezania, K. Shankar, and J. A. Tuszynski, Frontiers in Molecular Biosciences 8, 150 (2021b).
  * Schnyder  _et al._ (2009) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, AIP Conference Proceedings 1134, 10 (2009), https://aip.scitation.org/doi/pdf/10.1063/1.3149481 .
  * Delplace  _et al._ (2011) P. Delplace, D. Ullmo, and G. Montambaux, Phys. Rev. B 84, 195452 (2011).
  * Zhu  _et al._ (2019) L. Zhu, E. Prodan, and K. H. Ahn, Phys. Rev. B 99, 041117 (2019).
  * Li  _et al._ (2017) C. Li, S. Lin, G. Zhang, and Z. Song, Physical Review B 96 (2017).
  * Hügel and Paredes (2014) D. Hügel and B. Paredes, Phys. Rev. A 89, 023619 (2014).
  * Padavić  _et al._ (2018) K. Padavić, S. S. Hegde, W. DeGottardi, and S. Vishveshwara, Phys. Rev. B 98, 024205 (2018).
  * Lieu  _et al._ (2020) S. Lieu, M. McGinley, and N. R. Cooper, Phys. Rev. Lett. 124, 040401 (2020).
  * McGinley and Cooper (2019) M. McGinley and N. R. Cooper, Phys. Rev. Research 1, 033204 (2019).
  * Campos Venuti  _et al._ (2017) L. Campos Venuti, Z. Ma, H. Saleur, and S. Haas, Phys. Rev. A 96, 053858 (2017).
  * Pérez-González  _et al._ (2019) B. Pérez-González, M. Bello, A. Gómez-León, and G. Platero, Phys. Rev. B 99, 035146 (2019).
  * Vemu  _et al._ (2017) A. Vemu, J. Atherton, J. O. Spector, C. A. Moores, and A. Roll-Mecak, Mol Biol Cell 28, 3564 (2017).
  * Prodan (2011) E. Prodan, Journal of Physics A: Mathematical and Theoretical 44, 113001 (2011).
  * Claes and Hughes (2020) J. Claes and T. L. Hughes, Phys. Rev. B 101, 224201 (2020).
  * Mondragon-Shem  _et al._ (2014) I. Mondragon-Shem, T. L. Hughes, J. Song, and E. Prodan, Phys. Rev. Lett. 113, 046802 (2014).
  * Gadadhar  _et al._ (2017) S. Gadadhar, S. Bodakuntla, K. Natarajan, and C. Janke, Journal of Cell Science 130, 1347 (2017), https://journals.biologists.com/jcs/article-pdf/130/8/1347/1946365/jcs199471.pdf .