Neonatal seizures are sudden events in brain activity with detrimental effects in neurological functions usually related to epileptic fits. Though neonatal seizures can be identified from electroencephalography (EEG), this is a challenging endeavour since expert visual inspection of EEG recordings is time consuming and prone to errors due the data's nonstationarity and low signal-to-noise ratio. Towards the greater aim of automatic clinical decision making and monitoring, we propose a multi-output Gaussian process (MOGP) framework for neonatal EEG modelling. In particular, our work builds on the multi-output spectral mixture (MOSM) covariance kernel and shows that MOSM outperforms other commonly-used covariance functions in the literature when it comes to data imputation and hyperparameter-based seizure detection. To the best of our knowledge, our work is the first attempt at modelling and classifying neonatal EEG using MOGPs. Our main contributions are: i) the development of an MOGP-based framework for neonatal EEG analysis; ii) the experimental validation of the MOSM covariance kernel on real-world neonatal EEG for data imputation; and iii) the design of features for EEG based on MOSM hyperparameters and their validation for seizure detection (classification) in a patient specific approach.