we estimate the function $g_i$ by decomposing it as a linear combination of basis functions, which is chosen according to the application. The parameters of the basis functions are obtained by performing a $10$-fold cross-validation with $90\%$ training and $10\%$ test [@shandilya2011; @james2013] (see Methods for details). The particular technique used in this step is not essential as the dynamics will be low-dimensional. Other techniques such as compressive sensing based techniques [@brunton2015; @mangan2016; @wang2011] could be also employed.

-   [**Isolated dynamics and effective coupling.**]{} Since the network is heterogeneous, it has many low degree nodes for which $g_i$ is close to the uncoupled dynamics $f_i\approx f$. We first run a model-free estimation to coarsely classify nodes into low-degree nodes and hubs and estimate $f$ (see Supplementary Material). We then use the expression of $g_i$ recovered at low degree nodes to obtain an approximation for $f$, while $g_i$ recovered at hub nodes allows to estimate $\beta_i v$. We estimate the parameter $\beta_i$ by using Bayesian inference (see Methods).

-   [**Network structure and communities.**]{} Since $\beta_i=\alpha k_i$, we can recover the network’s degree distribution from knowledge of $\beta_i$ at every node. To reconstruct the community structures, we use the fact that the 