Neuromodulation with low-intensity focused ultrasound (LIFUS) holds significant promise for noninvasive treatment of neurological disorders, but its success relies heavily on accurately targeting specific brain regions. Computational model predictions can be used to optimize LIFUS, but uncertain acoustic tissue properties can affect prediction accuracy. The Monte Carlo method is often used to quantify the impact of uncertainties, but many iterations are generally needed for accurate estimates. We studied a surrogate model based on polynomial chaos expansion (PCE) to quantify the uncertainty in the LIFUS acoustic intensity field caused by tissue acoustic property uncertainties. The PCE approach was benchmarked against Monte Carlo method for LIFUS in three different head models. We also investigated the effect of the number of PCE samples on the accuracy of the surrogate model. Our results show that the PCE surrogate model requires only 20 simulation samples to estimate the mean and standard deviation of the acoustic intensity field with high accuracy compared to 100 samples needed for Monte Carlo method. The root mean squared percentage error (RMSPE) in the mean acoustic intensity field was less than 1.5%, with a maximum error of less than 0.5 W/cm