Feed-forward neural networks (NNs) are a staple machine learning method widely used in many areas of science and technology, including physical chemistry, computational chemistry, and materials informatics. While even a single-hidden-layer NN is a universal approximator, its expressive power is limited by the use of simple neuron activation functions (such as sigmoid functions) that are typically the same for all neurons. More flexible neuron activation functions would allow the use of fewer neurons and layers and thereby save computational cost and improve expressive power. We show that additive Gaussian process regression (GPR) can be used to construct optimal neuron activation functions that are individual to each neuron. An approach is also introduced that avoids nonlinear fitting of neural network parameters by defining them with rules. The resulting method combines the advantage of robustness of a linear regression with the higher expressive power of an NN. We demonstrate the approach by fitting the potential energy surfaces of the water molecule and formaldehyde. Without requiring any nonlinear optimization, the additive-GPR-based approach outperforms a conventional NN in the high-accuracy regime, where a conventional NN suffers more from overfitting.