We adopt a volume penalizing immersed boundary method for the simulation of pulsatile blood flow inside cerebral aneurysms. We show that the flow undergoes a transition from an orderly state at low physiological Reynolds numbers, in which the pulsatile forcing is closely followed in time, to a complex response with strongly increased high-frequency components at higher physiological Reynolds numbers, i.e., at higher flow rates and larger aneurysm sizes. The flow is computed by solving the Navier–Stokes equations for incompressible flow. Geometric complexity of aneurysms in the cerebrovascular system is captured by defining the fluid and solid domains using a so-called binary ‘masking function’, which is a key element in the immersed boundary method. The pulsatile variation of the flow rate is represented in terms of measured cross-sectionally averaged velocities in the vicinity of the aneurysm, obtained by noninvasive Transcranial Doppler sonography. Transition of the flow is found to arise in qualitatively the same way at all locations near the aneurysm bulge, quite independent of the solution component that is monitored. The numerical reliability of the predicted transition is quantified on the basis of practical upper and lower bounding solutions, expressing the sensitivity of the flow to uncertainties in the aneurysm geometry. We compute the spectrum of the response of the flow at various locations and flow conditions and quantify the transition in local pressure and velocity. The significant increase of small-scale, high-frequency structures at higher Reynolds numbers may have potential for clinical screening application in the future.