We propose a mathematical model based on multivariate Hawkes processes for the activity of (natural) neural networks, where the excitation function of each neuron is modelled as a 2nd order causal exponential spline. Then, we estimate the parameters of the multivariate Hawkes process given a realization using maximum likelihood estimation, which is numerically implemented using an adaptation of the Polyatomic Frank-Wolfe algorithm, selected for its sparsity-promoting properties. Finally, we test our algorithm on real spike train data, which consists of intracellular and extracellular recordings of 250 neurons in turtle brains.