Filter Implementation
There are a few different ways to proceed. They involve doing some algebra but it is more an issue of tedium. The notation of your analog filters is a little non standard, essentially your $p$ is $s$ in most books. 

Traditionally, phase wasn't considered important in hearing so standards tended to be specified by magnitude, so there is some ambiguity with respect to the digital implementation's phase response. You have 2 basic simple choices, a digital IIR filter or FIR filter. FIR is likely to be easier but easy and better are typically not the same. FIR can be direct or FFT based.  IIR has a number of choices such as biquads or direct. Floating point math has its own easy/better issues. 

One path way is to focus on the terms inside the magnitude brackets, i.e. the complex analog transfer functions. The other path is to focus on the analog magnitude.  

So, given the analog transfer functions, you can either use a table of Laplace transforms and deduce the impulse response of each component you listed, or use the bilinear transform to express each component that is in $s$ in terms of $z$. 

If you derive the impulse responses, you can sample points of the impulse response and provided a generally exponential decay, use those samples as coefficients of a FIR filter of sufficient length. The FIR filters are then applied sequentially. One nice property of ideal LTI filters, the order that you apply the filters are ideally interchangeable.  

If you choose the bilinear transform and do the algebra to derive the discrete time numerators and denominators, you need to be aware that the bilinear transform squeezes the response as it approaches the Nyquist frequency. One needs to plot the discrete time response to asses how the corresponding analog response is distorted. If some modest tweaking results in an acceptable response then, the numerators and denominators need to be factored so that they are compatible with the IIR implementation, such as a cascade of biquads.

If you go with the magnitudes, the task becomes a direct sampling from the analog frequency domain to the discrete time periodic frequency domain. One can choose a linear phase response, and appropriately impose symmetry and define the frequency domain FIR filter. This might not sound like an analog implementation.

The essential challenge is mapping a frequency response in $-\infty &lt; \Omega &lt; \infty$ to $-\pi \le \omega \le \pi$ 

There are other approaches like state variables, so this is not the entire story.