{'fraction_non_alphanumeric': 0.03375254928619986, 'fraction_numerical': 0.012372535690006799, 'mean_word_length': 4.6918165989553104, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Graph Prolongation Convolutional Networks: Explicitly Multiscale Machine Learning on Graphs, with Applications to Modeling of Biological Systems

We define a novel type of ensemble Graph Convolutional Network (GCN) model. Using optimized linear projection operators to map between spatial scales of graph, this ensemble model learns to aggregate information from each scale for its final prediction. We calculate these linear projection operators as the infima of an objective function relating the structure matrices used for each GCN. Equipped with these projections, our model (a Graph Prolongation-Convolutional Network) outperforms other GCN ensemble models at predicting the potential energy of monomer subunits in a coarse-grained mechanochemical simulation of microtubule bending. We demonstrate these performance gains by measuring an estimate of the FLOPs spent to train each model, as well as wall-clock time. Because our model learns at multiple scales, it is possible to train at each scale according to a predetermined schedule of coarse vs. fine training. We examine several such schedules adapted from the Algebraic Multigrid (AMG) literature, and quantify the computational benefit of each. Finally, we demonstrate how under certain assumptions, our graph prolongation layers may be decomposed into a matrix outer product of smaller GCN operations.

We define a novel type of ensemble Graph Convolutional Network (GCN) model. Using optimized linear projection operators to map between spatial scales of graph, this ensemble model learns to aggregate information from each scale for its final prediction. We calculate these linear projection operators as the infima of an objective function relating the structure matrices used for each GCN. Equipped with these projections, our model (a Graph Prolongation-Convolutional Network) outperforms other GCN ensemble models at predicting the potential energy of monomer subunits in a coarse-grained mechanochemical simulation of microtubule bending. We demonstrate these performance gains by measuring an estimate of the FLOPs spent to train each model, as well as wall-clock time. Because our model learns at multiple scales, it is possible to train at each scale according to a predetermined schedule of coarse vs. fine training. We examine several such schedules adapted from the Algebraic Multigrid (AMG) literature, and quantify the computational benefit of each. Finally, we demonstrate how under certain assumptions, our graph prolongation layers may be decomposed into a matrix outer product of smaller GCN operations.

Convolution and Graph Convolution
Recent successes of deep learning have demonstrated that the inductive bias of Convolutional Neural Networks (CNNs) makes them extremely efficient for analyzing data with an inherent grid structure, such as images or video.
In particular, many applications use these models to make per-node (per-pixel) predictions over grid graphs: examples include image segmentation, optical flow prediction, anticipating motion of objects in a scene, and facial detection/identification. Further work applies these methods to emulate physical models, by discretizing the input domain. Computational Fluid Dynamics and other scientific tasks featuring PDEs or ODEs on a domain discretized by a rectangular lattice have seen recent breakthroughs applying machine learning models, like CNNs to handle data which is structured this way. These models learn a set of local filters whose size is much smaller than the size of the domainthese filters may then be applied simultaneously across the entire domain, leveraging the fact that at a given scale the local behavior of the neighborhood around a pixel (voxel) is likely to be similar at all grid points.
Graph Convolutional Networks (GCNs) are a natural extension of the above idea of image 'filters' to arbitrary graphs rather than nD grids, which may be more suitable in some scientific contexts. Intuitively, GCNs replace the image filtering operation of CNNs with repeated passes of: 1) aggregation of information between nodes according to some structure matrix 2) nonlinear processing of data at each node according to some rule (most commonly a flat neural network which takes as separate input(s) the current vector at each node). We refer the reader to a recent survey by Bacciu et al (2019) for a more complete exploration of the taxonomy graph neural networks.

Microtubules
As an example of a dataset whose underlying graph is not a grid, we consider a coarse-grained simulation of a microtubule. Microtubules (MTs) are self-assembling nanostructures which are ubiquitous in living cells. MTs play important structural roles during cell division, cell growth, and separation of chromosomes (in eukaryotic cells) (Chakrabortty et al., 2018). Microtubules are comprised of a lattice structure of two conformations (α and β) of tubulin. Free-floating tubulin monomers associate energetically into dimer subunits, which then associate head-to-tail to form long chainlike complexes called protofilaments. Protofilaments as-sociate side-to side in a sheet; at some critical number of protofilaments (which varies between species and cell type) the sheet wraps closed to form a repeating helical lattice with a seam. See (Pampaloni & Florin, 2008), Page 303, Figure 1. Key properties of microtubules are: Dynamic instability: microtubules grow from one end by attracting free-floating tubulin monomers (VanBuren et al., 2005). Microtubules can spontaneously enter a "catastrophe" phase, in which they rapidly unravel, but can also "rescue" themselves from the catastrophe state and resume growth (Gardner et al., 2013;Shaw et al., 2003). Interactions: Microtubules interact with one another: they can dynamically avoid one another during the growth phase, or collide and bundle up, or collide and enter catastrophe (Tindemans et al., 2014). The exact mechanism governing these interactions is an area of current research. Structural strength: microtubules are very stiff, with a Young's Modulus estimated at ≈1GPa for some cases (Pampaloni & Florin, 2008). This stiffness is thought to play a role in reinforcing cell walls (Kis et al., 2002).
In this work we introduce a model which learns to reproduce the dynamics of a graph signal (defined as an association of each node in the network with a vector of discrete or realvalued labels) at multiple scales of graph resolution. We apply this model framework to predict the potential energy of each tubulin monomer in a mechanochemical simulation of a microtubule.

Simulation of MTs and Prior Work
Non-continuum, non-event-based simulation of large molecules is typically done by representing some molecular subunit as a particle/rigid body, and then defining rules for how these subunits interact energetically. Molecular Dynamics (MD) simulation is an expansive area of study and a detailed overview is beyond the scope of this paper. We instead proceed to describe in general terms some basic ideas relevant to the numerical simulation detailed in Section 3.1. MD simulations proceed from initial conditions by computing the forces acting on each particle (according to the potential energy interactions and any external forces, as required), determining their instantaneous velocities and acceleration accordingly, and then moving each particle by the distance it would move (given its velocity) for some small timestep. Many variations of this basic idea exist. The software we use for our MD simulations, LAMMPS (Plimpton, 1993) allows for many different types of update step: we use Verlet integration (updating particle position according to the central difference approximation of acceleration (Verlet, 1967)) and Langevin dynamics (modeling the behavior of a viscous surrounding solvent implicitly (Schneider & Stoll, 1978)). We also elect to use the microcanonical ensemble (NVE) -meaning that the update steps of the system maintain the total number of particles, the volume of the system, and the total energy (kinetic + potential). For more details of our simulation, see Section 3.1 and the source code, available in the Supplementary Material accompanying this paper. Independent of implementation details, a common component of many experiments in computational molecular dynamics is the prediction of the potential energy associated with a particular conformation of some molecular structure. Understanding the energetic behavior of a complex molecule yields insights into its macro-scale behavior: for instance, the problem of protein folding can be understood as seeking a lower-energy configuration. In this work, we apply graph convolutional networks, trained via a method we introduce, to predict these energy values for a section of microtubule.

Mathematical Background and Notation
Definitions: For all basic terms (graph, edge, vertex, degree) we use standard definitions. We use the notation {x i } b i=a to represent the sequence of x i indexed by the integers a, a + 1, a + 2, . . . b. Graph Laplacian: The graph Laplacian is the matrix given by is the adjacency matrix of G, and 1 is an appropriately sized vector of 1s. The graph Laplacian is given by some authors as the opposite sign. Graph Diffusion Distance (GDD): Given two graphs G 1 and G 2 , with |G 1 | ≤ |G 2 | the Graph Diffusion Distance D(G 1 , G 2 ) is given by: where C(P ) represents some set of constraints on P , and || · || F represents the Frobenius norm. We take C(P ) to be orthogonality: P T P = I. Note that since in general P is a rectangular matrix, it may not be the case that P P T = I. Our recent work (Author & Author, 2019) has examined variants of this distance measure, and techniques for efficiently calculating this distance. Detailing these is outside of the scope of this paper; all P matrices detailed in this work were calculated using the constrained optimization package Pymanopt (Townsend et al., 2016). Prolongation matrix: we use the term "prolongation matrix" to refer to a matrix which is the optimum of the minimization given in the definition of the GDD.

Model Architecture
The model we propose is an ensemble of GCNs at multiple scales, with optimized projection matrices performing the mapping in between scales (i.e. between ensemble members). More formally, be their structure matrices (for some chosen method z of calculating the structure matrix given the graph). In all experiments in this paper, we take z(G) = L(G), the graph Laplacian, as previously defined 1 . In an ensemble of Graph Convolutional Networks, let represent the parameters (filter matrix and bias vector) in layer l of the ith network.
We follow the GCN formulation given by Kipf and Welling (2016). Assuming an input tensor X of dimensions n × F (where n is the number of nodes in the graph and F is the dimension of the label at each node), we inductively define the layerwise update rules for a graph convolutional network as: where g m is the activation function of the mth layer.
When i = j − 1, let P i,j be an optimal (in either the sense of Graph Diffusion Distance, or in the sense we detail in section 3.
Our multiscale ensemble model is then constructed as: This model architecture is illustrated in Figure 1. When the P matrices are constant/fixed, we will refer to this model as a GPCN, for Graph Prolongation-Convolutional Network. However, we find in our experiments in Section 3.3 that validation error is further reduced when the P operators are tuned during the same gradient update step which updates the filter weights, which we refer to as an "adaptive" GPCN or A-GPCN. We explain our method for choosing Z i and optimizing P matrices in Section 3.3.

Dataset
In this Section we detail the process for generating the simulated microtubule data for comparison of our model with other GCN ensemble models. Our microtubule structure has 13 protofilaments (each 48 tubulin monomers long).
As in a biological microtubule, each tubulin monomer is Figure 1. Schematic of GPCN model. Data matrix X is fed into the model and repeatedly coarsened using optimized projection matrices P ik . These coarsened data matrices are separately fed into GCN models. The final output of the ensemble is the projected sum of the outputs of each component GCN.
offset (along the axis parallel to the protofilaments) from its neighbors in adjacent protofilaments, resulting in a helical structrure with a pitch of 3 tubulin units. We refer to this pitch as the "offset" in Section 3.2. Each monomer subunit (624 total) is represented as a point mass of 50 Dalton (8.30 × 10 −15 ng). The diameter of the whole structure is 26nm, and the length is ≈ 260nm. The model itself was constructed using Moltemplate (Jewett et al., 2013), a tool for constructing large regular molecules to be used in LAMMPS simulations.
For this model, we define energetic interactions for angles and associations only. No steric or dihedral interactions were used: for dihedrals, this was because the lattice structure of the tube meant any set of four molecules contributed to multiple, contradictory dihedral interactions 2 . Interaction energy of an association b was calculated using the "harmonic" bond style in LAMMPS, i.e.
where b 0 is the resting length and k is the strength of that interaction. The energy of an angle φ was similarly calculated using the "harmonic" angle style, where φ 0 is the resting angle and k is again the interaction strength. The resting lengths and angles for all energetic interactions were calculated using the resting geometry of our microtubule graph G mt : a LAMMPS script was used to print the value of every angle interaction in the model, and these were collected and grouped based on value (all 153 • angles, all 102• angles, etc). Each strength parameter was varied over the values in {.1, .3, .6, 1.0, 1.3, 1.6, 1.9, 2.0}, producing 8 5 parameter combinations. Langevin dynamics were used, but with small temperature, to ensure stability and emphasize mechanical interactions. See Table 1 and Figure 3 for details on each strength parameter.
GNU Parallel (Tange, 2011) was used to run a simulation for each combination of interaction parameters, using the particle dynamics simulation engine LAMMPS. In each simulation, we clamp the first two rings of tubulin monomers (nodes 1-26) in place, and apply force (in the negative y direction) to the final two rings of monomers (nodes 599-624). This force starts at 0 and ramps up during the first 128000 timesteps (one step = .5ns) to its maximum value of 3 × 10 −15 N. Once maximum force is reached, the simulation runs for 256000 additional timesteps, which in our experience was long enough for all particles to come to rest. See Figure 2 for an illustration (visualized with Ovito (Stukowski, 2010)) of the potential energy per-particle at the final frame of a typical simulation run. Every K = 32000 timesteps, we save the following for every particle: the position x, y, z; components of velocity v x , v y , v z ; components of force F x , F y , F z ; and the potential energy of the particle E. The dataset is then a concatenation of the 12 saved frames from every simulation run, comprising all combinations of input parameter values, where for each frame we have: x i , the input graph signal, a 624 × 13 matrix holding the position, velocity, and force on each particle, as well as values of the five interaction coefficients; and y i , the output graph signal, a 624 × 1 matrix holding the potential energy calculated for each particle.
During training, after a training/validation split, we normalize the data by taking the mean and standard deviation of the N train × 624 × 13 input and N train × 624 × 1 output tensors along their first axis. Each data tensor is then reduced by the mean and divided by the standard deviation so that all 624 × 13 inputs to the network have zero mean and unit standard deviation. We normalize using the training data only.

Graph Coarsening
In this Section we outline a procedure for determining the coarsened structure matrices to use in the hierarchy of GCN models comprising a GPCN. We use our microtubule graph  Table 1 (using the particle labels in this image) is applied everywhere in the model where that arrangement of particle and association types occurs in that position.
as an example. In this case, we have two a-priori guidelines for producing the reduced-order graphs: 1) the reduced models should still be a tube and 2) it makes sense from a biological point of view to coarsen by combining the α-β pairs into single subunits. Given these restrictions, we can explore the space of coarsened graphs and find the coarse graph which is nearest to our original graph (under the GDD).
Our microtubule model is a tube of length 48 units, 13 units per complete "turn", and with the seam offset by three units. We generalize this notion as follows: Let p be the offset, and k be the number of monomers in one turn of the tube, and n the number of turns of a tube graph G Tube(n,k,p) . The graph used in our simulation is thus G mt = G Tube(48,13,3) . We pick the medium scale model G inter to be G Tube(24,13,1) , as this is the result of combining each α-β pair of tubulin monomer units in the fine scale, into one tubulin dimer unit in the medium scale. We pick the coarsest graph G coarse by searching over possible offset tube graphs. Namely, we vary k ∈ {3, 4, . . . 12} and p ∈ {0, 1, 2, 3}, and compute the optimal P * and its associated distance D(G Tube(24,k,p) , G mt |P = P * ). Figure 4 shows the distance between G mt and various other tube graphs as parameters p and k are varied. The nearest G Tube(24,k,p) to G mt is that with p = 0 and k = 3. Note that Figure 4 has two columns for each value of k: these represent the coarse edges along the seam having weight (relative to the other edges) 1 (marked with an S) or having weight 2 (no S). This is motivated by the fact that our initial condensing of each dimer pair condensed pairs of seam edges into single edges.

Comparison to Other GCN Ensemble Models
In this experiment we demonstrate the efficiency advantages of our model by comparing our approach to other ensemble Graph Convolutional Networks. Within each ensemble, each GCN model consists of several graph convolution lay-  Figure 4. Directed Graph Diffusion Distance (GDD) between offset tube graphs and Gmt. Table cells colored by value. We see from this comparison that the two graphs which are closest to Gmt are G Tube(24,3,0) and G Tube(24,3,0) with an edge weight of 2 for connections along the seam, motivating our choice of G Tube(24,3,0) (unweighted) as the coarsest graph in our hierarchy. ers, followed by several dense layers which are applied to each node separately (node-wise dense layers can be alternatively understood as a GCN layer with Z = I, although we implement it differently for efficiency reasons). The input to the dense layers is the node-wise concatenation of the output of each GCN layer. Each ensemble is the sum output of several such GCNs. We compare our models to 1, 2, and 3-member GCN ensembles with the same number of filters (but all using the original fine-scale structure matrix). For GPCN models, P matrices were calculated using Pymanopt (Townsend et al., 2016) to optimize Equation 1 subject to orthogonality constraints. The same P were used to initialize the (variable) P matrices of A-GPCN models. We also compare our model to the work of Abu-El-Haija et. al (2018), who introduce the N-GCN model: an ensemble GCN in which each ensemble member uses a different power Z r of the structure matrix (to aggregate information from neighborhoods of radius r). We include a N-GCN with radii (1,2,4) and a N-GCN with radii (1,2,4,8,16).
All models were trained with the same train/validation split, using ADAM with default hyperparameters, in TensorFlow (Abadi et al., 2016  eration in our models. The cost of a graph convolutional layer with n × n structure matrix Z, n × F input data X, and F × C filter matrix W is estimated as: nF (|Z| + C), where |Z| is the number of nonzero entries of Z. This is calculated as the sum of the costs of the two matrix multiplications X · W and Z · XW , with the latter assumed to be implemented as sparse matrix multiplication and therefore requiring O(|Z|nF ) operations. For implementation reasons, our GCN layers (across all models) do not use sparse multiplication; if support for arbitrary-dimensional sparse tensor outer products is included in TensorFlow in the future, we would expect the wall-clock times in Table 4 to decrease. The cost of a dense layer (with n × F input data X, and F × C filter matrix W ) applied to every node separately is estimated as: O(nF C). The cost of taking the dot product between a n × k matrix and a k × m matrix (for example, the restriction/prolongation by P ) is estimated as O(nmk).
We summarize the structure of each of our models in Table  2. In Figure 6 we show a comparison between each of these models, for one particular random seed (42). Error on the validation set is tracked as a function of computational cost expended to train the model (under our cost assumption given above). We see that all four GPCN models outperform the other types of ensemble model during early training, in the sense that they reach lower levels of error for the same amount of computational work performed. Additionally, Figure 6. Comparison of Normalized MSE on held-out validation data as a function of FLOPs expended for a variety of ensemble Graph Convolutional Network Models. We see that especially in early stages of training, our model formulation learns faster than an ensemble of 2, 3 or 5 GCNs with the same number of filters. The error plotted is the model's minimum error thus far (on the validation data). See supplementary material for a nonminimized version, and versions of same plots with wall-clock time substituted for estimated FLOP count.
the adaptive GPCN models outperform all other models in terms of absolute error: after the same number of training epochs (using the same random seed) they reach an order of magnitude lower error. Table 3 shows summary statistics for several runs of this experiment with varying random seeds; we see that the A-GPCN models consistently outperform all other models considered. Note that Figures 6,8, and 7 plot the Normalize Mean Squared Error (NMSE). This unitless value compares the output signal to the target after both are normalized by the procedure described in section 3.1.

Comparison: All-at-Once or Coarse-to-Fine Training
In this Section we compare the computational cost of training the entire GPCN at once, versus training the different 'resolutions' (meaning the different GCNs in the hierarchy) of the network according to a more complicated training schedule. This approach is motivated by recent work in coarse-to-fine training of both flat and convolutional neural networks (Scott & Mjolsness, 2019;Zhao et al., 2019;Haber et al., 2018;Dou & Wu, 2015;Ke et al., 2017), as well as the extensive literature on Algebraic MultiGrid (AMG) methods (Vaněk et al., 1996).
AMG solvers for differential equations on a mesh (which arises as the discretization of some volume to be simulated) proceed by performing numerical "smoothing steps" at multiple resolutions of discretization. The intuition behind this approach is that modes of error should be smooth at a spatial  scale which is equivalent to their wavelength, i.e. the solver shouldn't spend many cycles resolving long-wavelength errors at the finest scale, since they can be resolved more efficiently at the coarse scale. Given a solver and a hierarchy of discretizations, the AMG literature defines several types of training procedures or "cycle" types (F-cycle, V-cycle, Wcycle). These cycles can be understood as being specified by a recursion parameter γ, which controls how many times the smoothing or training algorithm visits all of coarser levels of the hierarchy in between smoothing steps at a given scale. For example, when γ = 1 the algorithm proceeds from fine to coarse and back again, performing one smoothing step at each resolution -a 'V' cycle.
We investigate the efficiency of training 3-level GPCN and A-GPCN (as described in Section 3.3), using multigrid-like training schedules with γ ∈ {0, 1, 2, 3}, as well as "coarse-to-fine" training: training the coarse model to convergence, then training the coarse and intermediate models together (until convergence), then finally training all three models at once. For coarse-to-fine training convergence was defined to have occurred once 10 epochs had passed without improvement of the validation error.
Our experiments (see Figure 7) show that these training schedules do result in a slight increase in efficiency of the GPCN model, especially during the early phase of training. However, we also find that these schedules are seemingly not compatible with the adaptive GPCN, as demonstrated by the fact that there are long periods of training with no improvement in validation loss. Notably, the coarse-to-fine schedule is an exception: the coarse-to-fine training of the A-GPCN outperformed even the GPCN with no multigrid training.

Comparison with DiffPool
Graph coarsening procedures are in general not differentiable. DiffPool (Ying et al., 2018) aims to address this by constructing an auxiliary GCN, whose output is a pooling matrix. Formally: Suppose that at layer l of a GCN we have a n l × n l structure matrix Z (l) and a n × F data matrix X (l) .
In addition to GCN layers as described in Section 2, Ying et. al define a pooling operation at layer l as: where GCN pool is an auxillary GCN with its own set of parameters θ , and σ is the softmax function. The output of GCN pool is a n × n coarse matrix, each row of which is softmaxed to produce an affinity matrix S whose rows each sum to 1, representing each fine-scale node being connected to one unit's worth of coarse-scale nodes. The coarsened structural and data matrices for the next layer are then calculated as: Clearly, the additional GCN layers required to produce S (l) incur additional computational cost. We compare our 3level GPCN (adaptive and not) models from the experiment in Section 3.3 to a model which has the same structure, but in which each P matrix is replaced by the appropriatelysized output of a DIFFPOOL module, and furthermore the coarsened structure matrices are produced as in Equation 2. We see that our GPCN model achieves comparable validation loss with less computational work, and our A-GPCN model additionally achieves lower absolute validation loss.

Differentiable Models of Molecular Dynamics
This work demonstrates the use of feed-forward neural networks to approximate the energetic potentials of a mechanochemical model of an organic molecule. Pertimestep, GCN models may not be as fast as highlyparallelized, optimized MD codes. However, neural networks are highly flexible function approximators: the GCN training approach outlined in this paper could also be used to train a GCN which predicts the energy levels per particle at the end of a simulation (once equilibrium is reached), given the boundary conditions and initial conditions of each particle. In the case of our MT experiments, approximately 3 × 10 5 steps were required to reach equilibrium. The computational work to generate a suitably large and diverse training set would then be amortized by the GCN's ability to generalize to initial conditions, boundary conditions, and hyperparameters outside of this data set. Furthermore, this GCN reduced model would be fully differentiable, making it possible to perform gradient descent with respect to any of these inputs.

Tensor Factorization
Recent work has re-examined GCNs in the context of the extensive literature on tensor decompositions. LanczosNet (Liao et al., 2019), uses QR decomposition of the structure matrix to aggregate information from large neighborhoods of the graph. The "Tensor Graph Convolutional Network" of Zhang et. al (2018), is a different decomposition method, based on graph factorization; a product of GCNs operating on each factor graph can be as accurate as a single GCN acting on the product graph. Since recent work (Scott & Mjolsness, 2019) has shown that the GDD of a graph product is bounded by the distances between the factor graphs, it seems reasonable to combine both ideas into a model which uses a separate GPCN for each factor. One major benefit of this approach would be that a transfer-learning style approach can be used. For example, we could train a product of two GCN models on a short section of microtubule; and then re-use the weights in a model that predicts energetic potentials for a longer microtubule. This would allow us to extend our approach to MT models whose lengths are biologically relevant, e.g. 10 3 tubulin monomers.

Conclusion
We introduce a new type of graph ensemble model which explicitly learns to approximate behavior at multiple levels of coarsening. Our model outperforms several other types of GCN, including both other ensemble models and a model which coarsens the original graph using DiffPool. We also explore the effect of various training schedules, discovering that A-GPCNs can be effectively trained using a coarse-to-fine training schedule. We present the first use of GCNs to approximate energetic potentials in a model of a microtubule.