There is no strict definition of suitability of an activation function for neural networks. Instead there are a number of desirable traits, and functions that don't meet them or come close enough may perform badly in general (but those functions may still work in specific cases)
If you are using gradient descent as a training method, then differentiability is closest to being a requirement. However, even then as you noticed with ReLU, it is not an absolute requirement, provided behaviour close to discontinuities is good. For example $\frac{1}{x}$ or $log{x}$ make bad choices here due to how they behave for values near to $0$.
Wikipedia summarises some important desirable traits, plus compares popular activation functions in this table.
Nonlinear – When the activation function is non-linear, then a two-layer neural network can be proven to be a universal function approximator. The identity activation function does not satisfy this property. When multiple layers use the identity activation function, the entire network is equivalent to a single-layer model.
Range – When the range of the activation function is finite, gradient-based training methods tend to be more stable, because pattern presentations significantly affect only limited weights. When the range is infinite, training is generally more efficient because pattern presentations significantly affect most of the weights. In the latter case, smaller learning rates are typically necessary.
Continuously differentiable – This property is desirable (RELU is not continuously differentiable and has some issues with gradient-based optimization, but it is still possible) for enabling gradient-based optimization methods. The binary step activation function is not differentiable at 0, and it differentiates to 0 for all other values, so gradient-based methods can make no progress with it.
Monotonic – When the activation function is monotonic, the error surface associated with a single-layer model is guaranteed to be convex.
Smooth functions with a monotonic derivative – These have been shown to generalize better in some cases.
Approximates identity near the origin – When activation functions have this property, the neural network will learn efficiently when its weights are initialized with small random values. When the activation function does not approximate identity near the origin, special care must be used when initializing the weights.
You can see that ReLU is an outlier on more than one of these traits. The reason that it is popular is that despite these weaknesses, the things it does well - speed of operation plus helping combat the vanishing gradient problem - make it a solid choice for very deep networks. In fact the success of ReLU has inspired use of a number of similar-looking activation functions that hope to keep its benefits whilst having more of these desirable traits.
Ultmately though, almost any non-linear function can be used successfully in a neural network. The more of the desirable traits it has, the more likely you can use it in general cases with existing approaches to training.
