A fluid-structure interaction model that couples viscoelastic fluid motion induced by collective behavior of cilia to detailed axoneme mechanics is used to investigate bounds on cilia structure and mucus properties that determine effective fluid clearance. Dynein forcing is represented by a stochastic walker model that responds to local ATP concentrations provided by a biochemical network model. Cilium axonemes are modeled by large-deflection finite elements representing microtubules and viscoelastic springs representing connecting elements (nexins, radial spokes). Cilium motion is coupled to a viscoelastic fluid computation that models gel-like behavior of the mucus as well as possible Newtonian behavior of the periciliary fluid layer. Behavior of the viscoelastic fluid is prescribed at a microscopic level to avoid using continuum viscoelastic models of questionable validity. A lattice-based technique based upon a variational formulation of the Fokker-Planck equation is used to describe the viscoelastic fluid dynamics. A lattice Boltzmann method is applied to capture forcing of the viscoelastic mucus layer by concurrent airflow. The overall model exhibits natural formation of metachronal waves due to phase coupling of the dynein motion. Adjoint density analysis and uncertainty quantification techniques are applied to assess the stability of the transport induced by metachronal waves to perturbations in dynein walker rates, axoneme element rigidity, and mucus gel-formation process. The goal is not only to assess the robustness of the metachronal transport process, but also to identify elements within the overall transport mechanism that are most promising targets for pharmaceutical treatment of ciliary dysfunction.