In this work, we present a nonlinear principal component analysis (PCA) that identifies underlying sources causing the expression of spatial modes or patterns of activity in neuroimaging time series where these sources can interact to produce second-order modes. This nonlinear PCA uses a neural network architecture that embodies a specific form for the mixing of sources that is based on a second-order approximation to any general nonlinear mixing. The modes obtained have a unique rotation and scaling that does not depend on the biologically implausible constraints adopted by conventional PCA. Interactions among sources render the expression of any mode or brain system sensitive to the expression of others. The example considers interactions among functionally specialized brain systems (using a fMRI study of colour and motion processing).