The main objective of this research work is to develop an effective mathematical model of cardiac conduction system using a heterogeneous whole-heart model. The model is in the form of a system of modified Van der Pol and FitzHugh-Nagumo differential equations capable of describing the heart dynamics. The proposed model extends the range of normal and pathological electrocardiogram (ECG) waveforms that can be generated by the model. The effects of the respiratory sinus arrhythmia (RSA) and the Mayer waves (MW) are both incorporated to modulate the intrinsic frequency of the main oscillator that represents the sinoatrial node. Also, three pathological conditions are incorporated into the model. The heart rate variability (HRV) phenomenon is incorporated into the synthetic ECGs produced which yields valuable information about the cardiovascular health and the performance of the autonomic nervous system. The spectral analysis of the generated RR tachogram delivers power spectrums that resemble those obtained from real recordings. Also, the proposed model generates synthetic ECGs that characteristic the three considered pathological conditions, namely, the tall T wave, the ECG with U wave, and the Wolff-Parkinson-White syndrome. In general, the significance of this research work is in developing a mathematical model that represents the interactions between different pacemakers and allows analysis of cardiac rhythms. To show the effectiveness and the accuracy of the presented model, the results are compared to published results. The proposed model can be a useful tool to study the influences of different physiological conditions on the profile of the ECG. The synthetic ECG signals produced can be used as signal sources for the assessment of diagnostic ECG signal processing devices.